Trait Anxiety and State Inference - Main analyses
Ondrej Zika
17/3/2023
Preparatory steps
- identify root folder
- restore computational environment from $ROOT/renv.lock
Load packages
required_packages = c("confintr", "tidyr","ggplot2", "effectsize", "easystats", "Jmisc", "sjmisc", "plyr", "lme4", "lmerTest", "emmeans", "dplyr", "ggpubr", "purrr", "broom", "plotrix", "VGAM", "reshape2", "PupillometryR")
#detach(package:lmerTest, unload=TRUE)
#library(scales)
#library(ggplot2)
lapply(required_packages, library, character.only = TRUE)##
## Attaching package: 'effectsize'## The following object is masked from 'package:confintr':
##
## oddsratio## # Attaching packages: easystats 0.6.0
## ✔ bayestestR 0.13.0 ✔ correlation 0.8.3
## ✔ datawizard 0.6.5 ✔ insight 0.19.0
## ✔ modelbased 0.8.6 ✔ performance 0.10.2
## ✔ parameters 0.20.2 ✔ report 0.5.6
## ✔ see 0.7.4##
## Attaching package: 'Jmisc'## The following object is masked from 'package:parameters':
##
## demean## The following object is masked from 'package:datawizard':
##
## demean## The following object is masked from 'package:ggplot2':
##
## %+%##
## Attaching package: 'sjmisc'## The following objects are masked from 'package:datawizard':
##
## center, empty_rows, remove_empty_rows, reshape_longer, to_factor,
## to_numeric## The following object is masked from 'package:tidyr':
##
## replace_na##
## Attaching package: 'plyr'## The following object is masked from 'package:here':
##
## here## Loading required package: Matrix##
## Attaching package: 'Matrix'## The following objects are masked from 'package:tidyr':
##
## expand, pack, unpack##
## Attaching package: 'lmerTest'## The following object is masked from 'package:lme4':
##
## lmer## The following object is masked from 'package:stats':
##
## step##
## Attaching package: 'dplyr'## The following objects are masked from 'package:plyr':
##
## arrange, count, desc, failwith, id, mutate, rename, summarise,
## summarize## The following object is masked from 'package:Jmisc':
##
## recode## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union##
## Attaching package: 'ggpubr'## The following object is masked from 'package:plyr':
##
## mutate##
## Attaching package: 'purrr'## The following object is masked from 'package:plyr':
##
## compact## The following object is masked from 'package:sjmisc':
##
## is_empty##
## Attaching package: 'plotrix'## The following object is masked from 'package:datawizard':
##
## rescale## Loading required package: stats4## Loading required package: splines##
## Attaching package: 'VGAM'## The following object is masked from 'package:tidyr':
##
## fill##
## Attaching package: 'reshape2'## The following object is masked from 'package:tidyr':
##
## smiths## Loading required package: rlang##
## Attaching package: 'rlang'## The following objects are masked from 'package:purrr':
##
## %@%, flatten, flatten_chr, flatten_dbl, flatten_int, flatten_lgl,
## flatten_raw, invoke, splice## The following object is masked from 'package:sjmisc':
##
## is_empty## [[1]]
## [1] "confintr" "knitr" "here" "pacman" "stats" "graphics"
## [7] "grDevices" "datasets" "utils" "methods" "base"
##
## [[2]]
## [1] "tidyr" "confintr" "knitr" "here" "pacman" "stats"
## [7] "graphics" "grDevices" "datasets" "utils" "methods" "base"
##
## [[3]]
## [1] "ggplot2" "tidyr" "confintr" "knitr" "here" "pacman"
## [7] "stats" "graphics" "grDevices" "datasets" "utils" "methods"
## [13] "base"
##
## [[4]]
## [1] "effectsize" "ggplot2" "tidyr" "confintr" "knitr"
## [6] "here" "pacman" "stats" "graphics" "grDevices"
## [11] "datasets" "utils" "methods" "base"
##
## [[5]]
## [1] "see" "report" "parameters" "performance" "modelbased"
## [6] "insight" "datawizard" "correlation" "bayestestR" "easystats"
## [11] "effectsize" "ggplot2" "tidyr" "confintr" "knitr"
## [16] "here" "pacman" "stats" "graphics" "grDevices"
## [21] "datasets" "utils" "methods" "base"
##
## [[6]]
## [1] "Jmisc" "see" "report" "parameters" "performance"
## [6] "modelbased" "insight" "datawizard" "correlation" "bayestestR"
## [11] "easystats" "effectsize" "ggplot2" "tidyr" "confintr"
## [16] "knitr" "here" "pacman" "stats" "graphics"
## [21] "grDevices" "datasets" "utils" "methods" "base"
##
## [[7]]
## [1] "sjmisc" "Jmisc" "see" "report" "parameters"
## [6] "performance" "modelbased" "insight" "datawizard" "correlation"
## [11] "bayestestR" "easystats" "effectsize" "ggplot2" "tidyr"
## [16] "confintr" "knitr" "here" "pacman" "stats"
## [21] "graphics" "grDevices" "datasets" "utils" "methods"
## [26] "base"
##
## [[8]]
## [1] "plyr" "sjmisc" "Jmisc" "see" "report"
## [6] "parameters" "performance" "modelbased" "insight" "datawizard"
## [11] "correlation" "bayestestR" "easystats" "effectsize" "ggplot2"
## [16] "tidyr" "confintr" "knitr" "here" "pacman"
## [21] "stats" "graphics" "grDevices" "datasets" "utils"
## [26] "methods" "base"
##
## [[9]]
## [1] "lme4" "Matrix" "plyr" "sjmisc" "Jmisc"
## [6] "see" "report" "parameters" "performance" "modelbased"
## [11] "insight" "datawizard" "correlation" "bayestestR" "easystats"
## [16] "effectsize" "ggplot2" "tidyr" "confintr" "knitr"
## [21] "here" "pacman" "stats" "graphics" "grDevices"
## [26] "datasets" "utils" "methods" "base"
##
## [[10]]
## [1] "lmerTest" "lme4" "Matrix" "plyr" "sjmisc"
## [6] "Jmisc" "see" "report" "parameters" "performance"
## [11] "modelbased" "insight" "datawizard" "correlation" "bayestestR"
## [16] "easystats" "effectsize" "ggplot2" "tidyr" "confintr"
## [21] "knitr" "here" "pacman" "stats" "graphics"
## [26] "grDevices" "datasets" "utils" "methods" "base"
##
## [[11]]
## [1] "emmeans" "lmerTest" "lme4" "Matrix" "plyr"
## [6] "sjmisc" "Jmisc" "see" "report" "parameters"
## [11] "performance" "modelbased" "insight" "datawizard" "correlation"
## [16] "bayestestR" "easystats" "effectsize" "ggplot2" "tidyr"
## [21] "confintr" "knitr" "here" "pacman" "stats"
## [26] "graphics" "grDevices" "datasets" "utils" "methods"
## [31] "base"
##
## [[12]]
## [1] "dplyr" "emmeans" "lmerTest" "lme4" "Matrix"
## [6] "plyr" "sjmisc" "Jmisc" "see" "report"
## [11] "parameters" "performance" "modelbased" "insight" "datawizard"
## [16] "correlation" "bayestestR" "easystats" "effectsize" "ggplot2"
## [21] "tidyr" "confintr" "knitr" "here" "pacman"
## [26] "stats" "graphics" "grDevices" "datasets" "utils"
## [31] "methods" "base"
##
## [[13]]
## [1] "ggpubr" "dplyr" "emmeans" "lmerTest" "lme4"
## [6] "Matrix" "plyr" "sjmisc" "Jmisc" "see"
## [11] "report" "parameters" "performance" "modelbased" "insight"
## [16] "datawizard" "correlation" "bayestestR" "easystats" "effectsize"
## [21] "ggplot2" "tidyr" "confintr" "knitr" "here"
## [26] "pacman" "stats" "graphics" "grDevices" "datasets"
## [31] "utils" "methods" "base"
##
## [[14]]
## [1] "purrr" "ggpubr" "dplyr" "emmeans" "lmerTest"
## [6] "lme4" "Matrix" "plyr" "sjmisc" "Jmisc"
## [11] "see" "report" "parameters" "performance" "modelbased"
## [16] "insight" "datawizard" "correlation" "bayestestR" "easystats"
## [21] "effectsize" "ggplot2" "tidyr" "confintr" "knitr"
## [26] "here" "pacman" "stats" "graphics" "grDevices"
## [31] "datasets" "utils" "methods" "base"
##
## [[15]]
## [1] "broom" "purrr" "ggpubr" "dplyr" "emmeans"
## [6] "lmerTest" "lme4" "Matrix" "plyr" "sjmisc"
## [11] "Jmisc" "see" "report" "parameters" "performance"
## [16] "modelbased" "insight" "datawizard" "correlation" "bayestestR"
## [21] "easystats" "effectsize" "ggplot2" "tidyr" "confintr"
## [26] "knitr" "here" "pacman" "stats" "graphics"
## [31] "grDevices" "datasets" "utils" "methods" "base"
##
## [[16]]
## [1] "plotrix" "broom" "purrr" "ggpubr" "dplyr"
## [6] "emmeans" "lmerTest" "lme4" "Matrix" "plyr"
## [11] "sjmisc" "Jmisc" "see" "report" "parameters"
## [16] "performance" "modelbased" "insight" "datawizard" "correlation"
## [21] "bayestestR" "easystats" "effectsize" "ggplot2" "tidyr"
## [26] "confintr" "knitr" "here" "pacman" "stats"
## [31] "graphics" "grDevices" "datasets" "utils" "methods"
## [36] "base"
##
## [[17]]
## [1] "VGAM" "splines" "stats4" "plotrix" "broom"
## [6] "purrr" "ggpubr" "dplyr" "emmeans" "lmerTest"
## [11] "lme4" "Matrix" "plyr" "sjmisc" "Jmisc"
## [16] "see" "report" "parameters" "performance" "modelbased"
## [21] "insight" "datawizard" "correlation" "bayestestR" "easystats"
## [26] "effectsize" "ggplot2" "tidyr" "confintr" "knitr"
## [31] "here" "pacman" "stats" "graphics" "grDevices"
## [36] "datasets" "utils" "methods" "base"
##
## [[18]]
## [1] "reshape2" "VGAM" "splines" "stats4" "plotrix"
## [6] "broom" "purrr" "ggpubr" "dplyr" "emmeans"
## [11] "lmerTest" "lme4" "Matrix" "plyr" "sjmisc"
## [16] "Jmisc" "see" "report" "parameters" "performance"
## [21] "modelbased" "insight" "datawizard" "correlation" "bayestestR"
## [26] "easystats" "effectsize" "ggplot2" "tidyr" "confintr"
## [31] "knitr" "here" "pacman" "stats" "graphics"
## [36] "grDevices" "datasets" "utils" "methods" "base"
##
## [[19]]
## [1] "PupillometryR" "rlang" "reshape2" "VGAM"
## [5] "splines" "stats4" "plotrix" "broom"
## [9] "purrr" "ggpubr" "dplyr" "emmeans"
## [13] "lmerTest" "lme4" "Matrix" "plyr"
## [17] "sjmisc" "Jmisc" "see" "report"
## [21] "parameters" "performance" "modelbased" "insight"
## [25] "datawizard" "correlation" "bayestestR" "easystats"
## [29] "effectsize" "ggplot2" "tidyr" "confintr"
## [33] "knitr" "here" "pacman" "stats"
## [37] "graphics" "grDevices" "datasets" "utils"
## [41] "methods" "base"#pacman::p_load(char = required_packages)
source(here::here("utils/r_functions.R"))Load additional tools
# load color palettes
pal <- get_colors("ond2")
pal2 <- get_colors("ond")
pal3 <- get_colors("ukr")
plot(NULL, xlim=c(0,length(pal2)), ylim=c(0,1),
xlab="", ylab="", xaxt="n", yaxt="n")
rect(0:(length(pal)-1), 0.7, 1:length(pal), 0.9, col=pal)
rect(0:(length(pal2)-1), 0.4, 1:length(pal2), 0.6, col=pal2)
rect(0:(length(pal3)-1), 0.1, 1:length(pal3), 0.3, col=pal3)
1. Analysis of ratings - stable cues
Stats
stable_data<-read.csv(here::here("data/stable_cues_data.csv"))
stable_data = assign_var_types(stable_data, c("trtype_str", "study_str", "id", "contingency", "ta"))
# estimate models
m1a<-lmer(prob ~ ta*contingency*trtype_str + (1|id) + (1|study_str), stable_data)## boundary (singular) fit: see help('isSingular')m1b<-lmer(prob ~ ta*contingency*trtype_str + (1|id) + (1+ta+contingency+trtype_str|study_str), stable_data)## boundary (singular) fit: see help('isSingular')# model comparison with and without random slope
anova(m1a, m1b)## refitting model(s) with ML (instead of REML)## Data: stable_data
## Models:
## m1a: prob ~ ta * contingency * trtype_str + (1 | id) + (1 | study_str)
## m1b: prob ~ ta * contingency * trtype_str + (1 | id) + (1 + ta + contingency + trtype_str | study_str)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## m1a 15 -232.54 -175.924 131.27 -262.54
## m1b 29 -200.53 -91.073 129.27 -258.54 0 14 1# winning model overview
anova(m1a)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## ta 0.0000 0.0000 1 309.818 0.0000 0.996861
## contingency 0.0078 0.0039 2 17.227 0.1455 0.865610
## trtype_str 11.7153 11.7153 1 308.261 435.8064 < 2.2e-16 ***
## ta:contingency 0.0270 0.0135 2 309.263 0.5023 0.605632
## ta:trtype_str 0.1857 0.1857 1 308.261 6.9064 0.009019 **
## contingency:trtype_str 1.8189 0.9094 2 308.261 33.8311 5.211e-14 ***
## ta:contingency:trtype_str 0.0781 0.0391 2 308.261 1.4528 0.235503
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# effectsize
effectsize::effectsize(anova(m1a), type="eta", alternative = "two.sided")## # Effect Size for ANOVA (Type III)
##
## Parameter | Eta2 (partial) | 95% CI
## ---------------------------------------------------------
## ta | 5.00e-08 | [0.00, 0.00]
## contingency | 0.02 | [0.00, 0.17]
## trtype_str | 0.59 | [0.52, 0.64]
## ta:contingency | 3.24e-03 | [0.00, 0.02]
## ta:trtype_str | 0.02 | [0.00, 0.06]
## contingency:trtype_str | 0.18 | [0.11, 0.25]
## ta:contingency:trtype_str | 9.34e-03 | [0.00, 0.04]# statistical assumptions
## homogeneity of variance
stable_data$residuals <- residuals(m1a)
stable_data$residuals_sq <- abs(residuals(m1a))^2
hist(stable_data$residuals_sq)
lev_m <- lm(residuals_sq ~ id, data=stable_data) #ANOVA of the squared residuals
anova(lev_m)## Analysis of Variance Table
##
## Response: residuals_sq
## Df Sum Sq Mean Sq F value Pr(>F)
## id 88 0.28219 0.0032067 1.267 0.08311 .
## Residuals 233 0.58972 0.0025310
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1levRes = car::leveneTest(residuals_sq ~ id, data=stable_data)
levRes## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 88 1.2283 0.1142
## 233## normality of residuals
ggplot(data=stable_data, aes(x=residuals)) +
geom_histogram()## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
car::qqPlot(stable_data$residuals)
## [1] 152 188# marginals and post-hocs
em1a1 = emmeans(m1a, specs = pairwise ~ trtype_str)## NOTE: Results may be misleading due to involvement in interactionsem1a1$emmeans## trtype_str emmean SE df lower.CL upper.CL
## hi-prob 0.718 0.023 1.23 0.528 0.909
## low-prob 0.302 0.023 1.23 0.112 0.493
##
## Results are averaged over the levels of: contingency
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95em1a1$contrasts## contrast estimate SE df t.ratio p.value
## (hi-prob) - (low-prob) 0.416 0.0199 242 20.876 <.0001
##
## Results are averaged over the levels of: contingency
## Degrees-of-freedom method: kenward-rogerem1a2 = emmeans(m1a, specs = pairwise ~ trtype_str*contingency)## NOTE: Results may be misleading due to involvement in interactionsem1a2$emmeans## trtype_str contingency emmean SE df lower.CL upper.CL
## hi-prob 60/40 0.608 0.0373 6.76 0.5195 0.697
## low-prob 60/40 0.417 0.0373 6.76 0.3285 0.506
## hi-prob 75/25 0.724 0.0185 6.24 0.6792 0.769
## low-prob 75/25 0.306 0.0185 6.24 0.2615 0.351
## hi-prob 90/10 0.822 0.0370 6.39 0.7327 0.911
## low-prob 90/10 0.184 0.0370 6.39 0.0945 0.273
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95em1a2$contrasts## contrast estimate SE df t.ratio p.value
## (hi-prob 60/40) - (low-prob 60/40) 0.1910 0.0387 242.1 4.932 <.0001
## (hi-prob 60/40) - (hi-prob 75/25) -0.1158 0.0393 49.2 -2.947 0.0522
## (hi-prob 60/40) - (low-prob 75/25) 0.3019 0.0393 49.2 7.686 <.0001
## (hi-prob 60/40) - (hi-prob 90/10) -0.2135 0.0384 243.7 -5.554 <.0001
## (hi-prob 60/40) - (low-prob 90/10) 0.4247 0.0384 243.7 11.047 <.0001
## (low-prob 60/40) - (hi-prob 75/25) -0.3068 0.0393 49.2 -7.810 <.0001
## (low-prob 60/40) - (low-prob 75/25) 0.1109 0.0393 49.2 2.823 0.0705
## (low-prob 60/40) - (hi-prob 90/10) -0.4046 0.0384 243.7 -10.523 <.0001
## (low-prob 60/40) - (low-prob 90/10) 0.2337 0.0384 243.7 6.078 <.0001
## (hi-prob 75/25) - (low-prob 75/25) 0.4177 0.0247 242.1 16.881 <.0001
## (hi-prob 75/25) - (hi-prob 90/10) -0.0977 0.0390 47.4 -2.507 0.1426
## (hi-prob 75/25) - (low-prob 90/10) 0.5405 0.0390 47.4 13.863 <.0001
## (low-prob 75/25) - (hi-prob 90/10) -0.5155 0.0390 47.4 -13.221 <.0001
## (low-prob 75/25) - (low-prob 90/10) 0.1228 0.0390 47.4 3.149 0.0318
## (hi-prob 90/10) - (low-prob 90/10) 0.6382 0.0381 242.1 16.732 <.0001
##
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 6 estimatesem1a3 = emmeans(m1a, specs = pairwise ~ trtype_str | contingency)## NOTE: Results may be misleading due to involvement in interactionsem1a3$emmeans## contingency = 60/40:
## trtype_str emmean SE df lower.CL upper.CL
## hi-prob 0.608 0.0373 6.76 0.5195 0.697
## low-prob 0.417 0.0373 6.76 0.3285 0.506
##
## contingency = 75/25:
## trtype_str emmean SE df lower.CL upper.CL
## hi-prob 0.724 0.0185 6.24 0.6792 0.769
## low-prob 0.306 0.0185 6.24 0.2615 0.351
##
## contingency = 90/10:
## trtype_str emmean SE df lower.CL upper.CL
## hi-prob 0.822 0.0370 6.39 0.7327 0.911
## low-prob 0.184 0.0370 6.39 0.0945 0.273
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95em1a3$contrasts## contingency = 60/40:
## contrast estimate SE df t.ratio p.value
## (hi-prob) - (low-prob) 0.191 0.0387 242 4.932 <.0001
##
## contingency = 75/25:
## contrast estimate SE df t.ratio p.value
## (hi-prob) - (low-prob) 0.418 0.0247 242 16.881 <.0001
##
## contingency = 90/10:
## contrast estimate SE df t.ratio p.value
## (hi-prob) - (low-prob) 0.638 0.0381 242 16.732 <.0001
##
## Degrees-of-freedom method: kenward-rogera = summary(em1a3$contrasts)
effectsize::t_to_eta2(t=a$t.ratio, df_error = a$df, alternative = "two.sided")## Eta2 (partial) | 95% CI
## -----------------------------
## 0.09 | [0.03, 0.17]
## 0.54 | [0.46, 0.61]
## 0.54 | [0.46, 0.60]# trend with anxiety
emtr<-emtrends(m1a, pairwise ~ trtype_str, var = "ta")## NOTE: Results may be misleading due to involvement in interactionssumtr<-summary(emtr$contrasts)
effectsize::t_to_eta2(t=sumtr$t.ratio, df_error = sumtr$df, type="eta", alternative = "two.sided")## Eta2 (partial) | 95% CI
## -----------------------------
## 0.03 | [0.00, 0.08]# one-way t-tests against true reinforcement (estimated by err)
df2 <- stable_data %>%
group_by(contingency, tabin, trtype_str) %>%
nest() %>%
dplyr::mutate(tt=map(data,~t.test(.x$err))) %>%
dplyr::mutate(tidied = map(tt, tidy)) %>%
unnest(tidied, .drop = T)
df2$p.value.fdr = p.adjust(df2$p.value, method = "fdr")
df2## # A tibble: 12 × 14
## # Groups: contingency, tabin, trtype_str [12]
## tabin contingency trtype_str data tt estimate statistic p.value parameter
## <fct> <fct> <fct> <lis> <lis> <dbl> <dbl> <dbl> <dbl>
## 1 lowA… 75/25 hi-prob <tib… <hte… -0.0535 -2.62 1.18e-2 47
## 2 lowA… 75/25 low-prob <tib… <hte… 0.0810 3.58 8.18e-4 47
## 3 hiAnx 75/25 hi-prob <tib… <hte… -0.0344 -1.32 1.95e-1 39
## 4 hiAnx 75/25 low-prob <tib… <hte… 0.0428 2.04 4.78e-2 39
## 5 hiAnx 60/40 hi-prob <tib… <hte… -0.0179 -0.424 6.77e-1 17
## 6 hiAnx 60/40 low-prob <tib… <hte… 0.0239 0.682 5.05e-1 17
## 7 hiAnx 90/10 hi-prob <tib… <hte… -0.0723 -2.07 5.42e-2 17
## 8 hiAnx 90/10 low-prob <tib… <hte… 0.0160 0.776 4.49e-1 17
## 9 lowA… 60/40 hi-prob <tib… <hte… -0.00191 -0.0421 9.67e-1 17
## 10 lowA… 60/40 low-prob <tib… <hte… 0.00247 0.0567 9.55e-1 17
## 11 lowA… 90/10 hi-prob <tib… <hte… -0.152 -3.51 2.52e-3 18
## 12 lowA… 90/10 low-prob <tib… <hte… 0.157 2.32 3.21e-2 18
## # … with 5 more variables: conf.low <dbl>, conf.high <dbl>, method <chr>,
## # alternative <chr>, p.value.fdr <dbl>effectsize::t_to_eta2(t=df2$statistic, df_error = df2$parameter, alternative = "two.sided")## Eta2 (partial) | 95% CI
## -----------------------------
## 0.13 | [0.01, 0.31]
## 0.21 | [0.04, 0.40]
## 0.04 | [0.00, 0.22]
## 0.10 | [0.00, 0.30]
## 0.01 | [0.00, 0.24]
## 0.03 | [0.00, 0.29]
## 0.20 | [0.00, 0.50]
## 0.03 | [0.00, 0.31]
## 1.04e-04 | [0.00, 0.03]
## 1.89e-04 | [0.00, 0.06]
## 0.41 | [0.07, 0.65]
## 0.23 | [0.00, 0.52]## Gather model predictions and summarize for plots below
stable_data$model_pred <- predict(m1a) Plots
Stables 1: High/Low cue
stable_data_df = stable_data %>%
group_by(id, trtype_str) %>%
summarise_at(c("model_pred","prob"), mean, na.rm = TRUE)
g <- ggplot(data = stable_data_df, aes(y = prob, x = trtype_str, fill = trtype_str)) +
geom_hline(yintercept = 0.75, color=pal[5], linetype="longdash") +
geom_hline(yintercept = 0.25, color=pal[7], linetype="longdash") +
geom_flat_violin(position = position_nudge(x = .2, y = 0), alpha = .8, lwd=1.2) +
geom_point(aes(y = prob, color=trtype_str), position = position_jitter(width = .15), size = 2, alpha = 0.5, show.legend = FALSE) +
geom_boxplot(width = .2, outlier.shape = NA, alpha = 0.7, lwd=1.2, show.legend = FALSE) +
stat_summary(geom="point", fun="mean", size=4, shape=23, stroke=1.5, aes(y=model_pred, group=trtype_str), color="black", show.legend = FALSE) +
expand_limits(x = 3) +
scale_color_manual(values = c(pal[5], pal[7])) +
scale_fill_manual(values = c(pal[5], pal[7]), name = "Stable cue", labels = c("High-prob Cue", "Low-prob Cue")) +
# coord_flip() +
theme_bw() +
raincloud_theme +
theme(legend.position = "right") +
labs(y= "Mean probability", x="Session") +
scale_x_discrete(labels=c("hi-prob" = "HC", "low-prob" = "LC")) +
guides(fill=guide_legend(nrow=2,byrow=TRUE))
g
ggsave(here::here("output/figures/Fig_stable_1.pdf"), width = 6, height = 5, dpi = 120)Stables 2: High/Low cue by session
stable_data_df2 = stable_data %>%
group_by(id, trtype_str, contingency) %>%
summarise_at(c("model_pred","prob"), mean, na.rm = TRUE)
w <- 0.4
g <- ggplot(data = stable_data_df2, aes(y = prob, x = contingency, fill = trtype_str))
h<-c(0.6, 0.75, 0.9)
#
for (s in seq(3)) {
g<- g + geom_segment(x=s-w,y=h[s], xend=s+w, yend=h[s] , linetype="dashed", color=pal[5], size=0.5 , alpha=0.3)
g <- g + geom_segment(x=s-w,y=1-h[s], xend=s+w, yend=1-h[s] , linetype="dashed", color=pal[7], size=0.5 , alpha=0.3)
}
g<- g+geom_flat_violin(position = position_nudge(x = .2, y = 0), alpha = .8, show.legend = FALSE, lwd=1.2) +
geom_point(aes(y = prob, colour=trtype_str), position = position_jitterdodge( jitter.width = .15, dodge.width = 0.2), size = 2, alpha = 0.8, show.legend = FALSE) +
geom_boxplot(width = .2, outlier.shape = NA, alpha = 0.7,lwd=1.2, show.legend = FALSE) +
stat_summary(geom="point", fun="mean", size=4, shape=23, aes(y=model_pred, group=trtype_str), color="black",
position = position_dodge( 0.2), stroke=1.5, show.legend = FALSE) +
expand_limits(x = 3) +
scale_color_manual(values = c(pal[5], pal[7])) +
scale_fill_manual(values = c(pal[5], pal[7])) +
theme_bw() +
raincloud_theme +
labs(y= "Mean probability", x="Session") +
guides(fill=guide_legend(nrow=2,byrow=TRUE))
g
ggsave(here::here("output/figures/Fig_stable_2.pdf"), width = 8, height = 5, dpi = 120)Stables 3: High/Low cue by median-split TA
stable_data_df3 = stable_data %>%
group_by(id, trtype_str, tabin) %>%
summarise_at(c("model_pred","prob"), mean, na.rm = TRUE)
g <- ggplot(data = stable_data_df3, aes(y = prob, x = trtype_str, fill = interaction(tabin,trtype_str))) +
geom_hline(yintercept = 0.75, color=pal[5], linetype="longdash") +
geom_hline(yintercept = 0.25, color=pal[7], linetype="longdash") +
geom_flat_violin(position = position_nudge(x = .2, y = 0), alpha = .8, lwd=1.2) +
geom_point(aes(color= interaction( tabin, trtype_str)), position = position_jitterdodge( jitter.width = .45, dodge.width = 0.05), size = 2, alpha = 0.5, show.legend=FALSE) +
geom_boxplot(width = .2, outlier.shape = NA, alpha = 0.7, lwd=1.2, show.legend = FALSE) +
stat_summary(geom="point", fun="mean", size=4, shape=23, aes(y=model_pred, group=interaction(tabin,trtype_str)), color="black",
position = position_dodge( 0.2), stroke=1.5, show.legend = FALSE) +
expand_limits(x = 3) +
scale_color_manual(values = c(pal[6], pal[5], pal[8], pal[7])) +
scale_fill_manual(values = c(pal[6], pal[5], pal[8], pal[7]), name = "", labels = c("HC: Low Anxiety", "HC: High Anxiety", "LC: Low Anxiety", "LC: High Anxiety")) +
theme_bw() +
raincloud_theme +
theme(legend.position = "right") +
labs(y= "Mean probability", x="Session") +
scale_x_discrete(labels=c("hi-prob" = "HC", "low-prob" = "LC")) +
guides(fill=guide_legend(nrow=4,byrow=TRUE))
g
ggsave(here::here("output/figures/Fig_stable_3.pdf"), width = 7, height = 5, dpi = 120)stable_data_df3 = stable_data %>%
group_by(id, trtype_str, ta_orig_scale) %>%
summarise_at(c("model_pred","prob"), mean, na.rm = TRUE)
g <- ggplot(data = stable_data_df3, aes(y = prob, x = ta_orig_scale, fill = trtype_str)) +
geom_hline(yintercept = 0.75, color=pal[5], linetype="longdash") +
geom_hline(yintercept = 0.25, color=pal[7], linetype="longdash") +
#geom_flat_violin(position = position_nudge(x = .2, y = 0), alpha = .8, lwd=1.2) +
geom_point(aes(color= trtype_str), position = position_jitterdodge( jitter.width = .45, dodge.width = 0.05), size = 2, alpha = 0.5, show.legend=FALSE) +
geom_smooth(aes(color= trtype_str, fill=trtype_str), method='lm', formula= y~x, alpha=0.3, show.legend = FALSE) +
#geom_boxplot(width = .2, outlier.shape = NA, alpha = 0.7, lwd=1.2, show.legend = FALSE) +
#stat_summary(geom="point", fun="mean", size=4, shape=23, aes(y=model_pred, group=interaction(tabin,trtype_str)), color="black",
# position = position_dodge( 0.2), stroke=1.5, show.legend = FALSE) +
expand_limits(x = 3) +
scale_color_manual(values = c(pal[5], pal[7], pal[8], pal[7])) +
scale_fill_manual(values = c(pal[5], pal[7], pal[8], pal[7]), name = "", labels = c("HC: Low Anxiety", "HC: High Anxiety", "LC: Low Anxiety", "LC: High Anxiety")) +
theme_bw() +
raincloud_theme +
theme(legend.position = "right") +
labs(y= "Mean probability", x="Trait anxiety") +
scale_x_discrete(labels=c("hi-prob" = "HC", "low-prob" = "LC")) +
guides(fill=guide_legend(nrow=4,byrow=TRUE)) +
xlim(20, 70)## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.g
ggsave(here::here("output/figures/Fig_stable_3_ver2.pdf"), width = 4, height = 5, dpi = 120)Stables 4: Error from true by cue type, session and TA
stable_data_df4 = stable_data %>%
group_by(id, trtype_str, contingency, tabin) %>%
summarise_at(c("model_pred", "err"), mean, na.rm = TRUE)
stable_data_df4 <- within(stable_data_df4, err[trtype_str=="hi-prob"] <- (err[trtype_str=="hi-prob"]))
# https://stackoverflow.com/questions/16026215/generate-ggplot2-boxplot-with-different-colours-for-multiple-groups
g <- ggplot(data = stable_data_df4, aes(y = err, x = trtype_str, fill = interaction( tabin, trtype_str))) +
geom_hline(yintercept=c(0), linetype="dashed") +
geom_point(aes(color= interaction( tabin, trtype_str)), position = position_jitterdodge( jitter.width = .45, dodge.width = 1),
size = 2, alpha = 0.5, show.legend=FALSE) +
geom_boxplot(inherit.aes = TRUE, width = 0.8, lwd=1.2, outlier.shape = NA, alpha = 0.7, show.legend = FALSE) +
stat_summary(geom="point", fun="mean", size=4, shape=23, aes(y=err, group=interaction(tabin,trtype_str)), color="black",
position = position_dodge( 0.8), stroke=1.5, show.legend = FALSE) +
expand_limits(x = 3) +
scale_color_manual(values = c(pal[6], pal[5], pal[8], pal[7])) +
scale_fill_manual(values = c(pal[6], pal[5], pal[8], pal[7]), name = "", labels = c("HC: Low Anxiety", "HC: High Anxiety", "LC: Low Anxiety", "LC: High Anxiety")) +
theme_bw() +
raincloud_theme +
theme(legend.position = "right",
strip.text.x = element_text(size = 12, face = "bold"),
strip.background = element_rect(color="black", fill="#ebebeb", size=0, linetype="solid")
) +
labs(y= "Error", x="Session") +
scale_x_discrete(labels=c("hi-prob" = "HC", "low-prob" = "LC")) +
ylim(-0.6, 0.6) +
facet_grid(cols=vars(contingency))
g
ggsave(here::here("output/figures/Fig_stable_4.pdf"), width = 8, height = 5, dpi = 120)2. Analysis of ratings - reversal cue
Stats
rev_data<-read.csv(here::here("data/reversal_cue_data.csv"))
rev_data = assign_var_types(rev_data, c("phase_str", "study_str", "id", "contingency", "ta"))
dt_count <- rev_data %>%
group_by(contingency) %>%
summarise(count = n_distinct(id))
dt_count## # A tibble: 3 × 2
## contingency count
## <fct> <int>
## 1 60/40 36
## 2 75/25 88
## 3 90/10 37m2a <-lmer(prob ~ ta*contingency*phase_str + (1|id) + (1|study_str), rev_data) ## boundary (singular) fit: see help('isSingular')m2b<-lmer(prob ~ ta*contingency*phase_str + (1|id) + (1+ta+contingency+phase_str|study_str), rev_data)## boundary (singular) fit: see help('isSingular')anova(m2a, m2b)## refitting model(s) with ML (instead of REML)## Data: rev_data
## Models:
## m2a: prob ~ ta * contingency * phase_str + (1 | id) + (1 | study_str)
## m2b: prob ~ ta * contingency * phase_str + (1 | id) + (1 + ta + contingency + phase_str | study_str)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## m2a 15 -156.7 -100.079 93.349 -186.7
## m2b 29 -128.7 -19.237 93.349 -186.7 0.0018 14 1# anova
anova(m2a)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## ta 0.0454 0.0454 1 86.976 1.6244 0.2058701
## contingency 0.0346 0.0173 2 276.994 0.6199 0.5387247
## phase_str 4.6362 4.6362 1 234.049 165.9236 < 2.2e-16 ***
## ta:contingency 0.0453 0.0226 2 277.980 0.8103 0.4457800
## ta:phase_str 0.3843 0.3843 1 234.049 13.7552 0.0002601 ***
## contingency:phase_str 0.4950 0.2475 2 234.049 8.8578 0.0001958 ***
## ta:contingency:phase_str 0.0922 0.0461 2 234.049 1.6495 0.1943649
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# statistical assumptions
## homogeneity of variance
rev_data$residuals <- residuals(m2a)
rev_data$residuals_sq <- abs(residuals(m2a))^2
hist(rev_data$residuals_sq)
lev_m <- lm(residuals_sq ~ id, data=rev_data) #ANOVA of the squared residuals
anova(lev_m)## Analysis of Variance Table
##
## Response: residuals_sq
## Df Sum Sq Mean Sq F value Pr(>F)
## id 88 0.08185 0.00093012 0.8141 0.8676
## Residuals 233 0.26620 0.00114247## normality of residuals
ggplot(data=rev_data, aes(x=residuals)) +
geom_histogram()## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
car::qqPlot(rev_data$residuals)
## [1] 284 151# effectsize
effectsize::effectsize(anova(m2a), type="eta", alternative = "two.sided")## # Effect Size for ANOVA (Type III)
##
## Parameter | Eta2 (partial) | 95% CI
## --------------------------------------------------------
## ta | 0.02 | [0.00, 0.11]
## contingency | 4.46e-03 | [0.00, 0.03]
## phase_str | 0.41 | [0.32, 0.50]
## ta:contingency | 5.80e-03 | [0.00, 0.03]
## ta:phase_str | 0.06 | [0.01, 0.12]
## contingency:phase_str | 0.07 | [0.02, 0.14]
## ta:contingency:phase_str | 0.01 | [0.00, 0.05]# post-hocs
em2a1 = emmeans(m2a, specs = pairwise ~ phase_str*contingency)## NOTE: Results may be misleading due to involvement in interactionsem2a1$emmeans## phase_str contingency emmean SE df lower.CL upper.CL
## acq 60/40 0.630 0.0386 7.42 0.539 0.720
## ext 60/40 0.483 0.0386 7.42 0.393 0.573
## acq 75/25 0.714 0.0209 4.64 0.659 0.769
## ext 75/25 0.455 0.0209 4.64 0.400 0.510
## acq 90/10 0.763 0.0382 7.02 0.673 0.853
## ext 90/10 0.384 0.0382 7.02 0.293 0.474
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95em2a1$contrasts## contrast estimate SE df t.ratio p.value
## (acq 60/40) - (ext 60/40) 0.1464 0.0395 228.6 3.707 0.0035
## (acq 60/40) - (acq 75/25) -0.0842 0.0394 80.1 -2.138 0.2790
## (acq 60/40) - (ext 75/25) 0.1743 0.0394 80.1 4.422 0.0004
## (acq 60/40) - (acq 90/10) -0.1337 0.0392 229.4 -3.407 0.0100
## (acq 60/40) - (ext 90/10) 0.2458 0.0392 229.4 6.265 <.0001
## (ext 60/40) - (acq 75/25) -0.2307 0.0394 80.1 -5.853 <.0001
## (ext 60/40) - (ext 75/25) 0.0278 0.0394 80.1 0.706 0.9807
## (ext 60/40) - (acq 90/10) -0.2801 0.0392 229.4 -7.139 <.0001
## (ext 60/40) - (ext 90/10) 0.0994 0.0392 229.4 2.534 0.1187
## (acq 75/25) - (ext 75/25) 0.2585 0.0252 228.6 10.247 <.0001
## (acq 75/25) - (acq 90/10) -0.0495 0.0390 77.6 -1.266 0.8021
## (acq 75/25) - (ext 90/10) 0.3301 0.0390 77.6 8.453 <.0001
## (ext 75/25) - (acq 90/10) -0.3079 0.0390 77.6 -7.886 <.0001
## (ext 75/25) - (ext 90/10) 0.0716 0.0390 77.6 1.833 0.4509
## (acq 90/10) - (ext 90/10) 0.3795 0.0389 228.6 9.759 <.0001
##
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 6 estimatesemstr <- summary(em2a1$contrasts)
effectsize::t_to_eta2(t=emstr$t.ratio, df_error = emstr$df, alternative = "two.sided")## Eta2 (partial) | 95% CI
## -----------------------------
## 0.06 | [0.01, 0.12]
## 0.05 | [0.00, 0.17]
## 0.20 | [0.06, 0.34]
## 0.05 | [0.01, 0.11]
## 0.15 | [0.07, 0.23]
## 0.30 | [0.14, 0.44]
## 6.19e-03 | [0.00, 0.08]
## 0.18 | [0.10, 0.27]
## 0.03 | [0.00, 0.08]
## 0.31 | [0.22, 0.40]
## 0.02 | [0.00, 0.12]
## 0.48 | [0.32, 0.60]
## 0.44 | [0.28, 0.57]
## 0.04 | [0.00, 0.16]
## 0.29 | [0.20, 0.38]# relationship with anxiety
jt<-joint_tests(m2a, by = "phase_str")
jt ## phase_str = acq:
## model term df1 df2 F.ratio p.value
## ta 1 154.55 0.904 0.3432
## contingency 2 112.09 5.375 0.0059
## ta:contingency 2 229.56 0.044 0.9573
##
## phase_str = ext:
## model term df1 df2 F.ratio p.value
## ta 1 154.55 9.426 0.0025
## contingency 2 112.09 4.760 0.0104
## ta:contingency 2 229.56 2.319 0.1006effectsize::F_to_eta2(f=jt$F.ratio, df=jt$df1, df_error = jt$df2, type="eta", alternative = "two.sided")## Eta2 (partial) | 95% CI
## -----------------------------
## 5.82e-03 | [0.00, 0.05]
## 0.06 | [0.01, 0.14]
## 0.09 | [0.01, 0.19]
## 0.08 | [0.00, 0.18]
## 3.83e-04 | [0.00, 0.00]
## 0.02 | [0.00, 0.06]emtr <- emtrends(m2a, pairwise ~ phase_str, var = "ta")## NOTE: Results may be misleading due to involvement in interactionsemtr## $emtrends
## phase_str ta.trend SE df lower.CL upper.CL
## acq 0.00160 0.00168 155 -0.00173 0.00493
## ext -0.00517 0.00168 155 -0.00850 -0.00184
##
## Results are averaged over the levels of: contingency
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## acq - ext 0.00677 0.00183 229 3.709 0.0003
##
## Results are averaged over the levels of: contingency
## Degrees-of-freedom method: kenward-rogersumtr<-summary(emtr$contrasts)
sumtr## contrast estimate SE df t.ratio p.value
## acq - ext 0.00677 0.00183 229 3.709 0.0003
##
## Results are averaged over the levels of: contingency
## Degrees-of-freedom method: kenward-rogereffectsize::t_to_eta2(t=sumtr$t.ratio, df_error = sumtr$df, type="eta", alternative = "two.sided")## Eta2 (partial) | 95% CI
## -----------------------------
## 0.06 | [0.01, 0.12]## there is no three-way interaction but this is for visualisation purpose
emtr <- emtrends(m2a, pairwise ~ phase_str * contingency, var = "ta")
emtr## $emtrends
## phase_str contingency ta.trend SE df lower.CL upper.CL
## acq 60/40 0.00146 0.00279 303 -0.00402 0.00694
## ext 60/40 -0.00499 0.00279 303 -0.01047 0.00049
## acq 75/25 0.00121 0.00196 264 -0.00264 0.00507
## ext 75/25 -0.00193 0.00196 264 -0.00578 0.00192
## acq 90/10 0.00213 0.00274 302 -0.00326 0.00751
## ext 90/10 -0.00859 0.00274 302 -0.01398 -0.00320
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## (acq 60/40) - (ext 60/40) 0.006454 0.00352 229 1.835 0.4456
## (acq 60/40) - (acq 75/25) 0.000247 0.00315 271 0.079 1.0000
## (acq 60/40) - (ext 75/25) 0.003391 0.00315 271 1.078 0.8899
## (acq 60/40) - (acq 90/10) -0.000664 0.00349 230 -0.190 1.0000
## (acq 60/40) - (ext 90/10) 0.010053 0.00349 230 2.880 0.0491
## (ext 60/40) - (acq 75/25) -0.006206 0.00315 271 -1.972 0.3610
## (ext 60/40) - (ext 75/25) -0.003062 0.00315 271 -0.973 0.9261
## (ext 60/40) - (acq 90/10) -0.007117 0.00349 230 -2.039 0.3236
## (ext 60/40) - (ext 90/10) 0.003600 0.00349 230 1.031 0.9070
## (acq 75/25) - (ext 75/25) 0.003144 0.00239 229 1.316 0.7757
## (acq 75/25) - (acq 90/10) -0.000911 0.00310 270 -0.294 0.9997
## (acq 75/25) - (ext 90/10) 0.009806 0.00310 270 3.159 0.0216
## (ext 75/25) - (acq 90/10) -0.004055 0.00310 270 -1.306 0.7814
## (ext 75/25) - (ext 90/10) 0.006662 0.00310 270 2.146 0.2668
## (acq 90/10) - (ext 90/10) 0.010717 0.00345 229 3.103 0.0259
##
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 6 estimatessumtr<-summary(emtr$contrasts)
sumtr## contrast estimate SE df t.ratio p.value
## (acq 60/40) - (ext 60/40) 0.006454 0.00352 229 1.835 0.4456
## (acq 60/40) - (acq 75/25) 0.000247 0.00315 271 0.079 1.0000
## (acq 60/40) - (ext 75/25) 0.003391 0.00315 271 1.078 0.8899
## (acq 60/40) - (acq 90/10) -0.000664 0.00349 230 -0.190 1.0000
## (acq 60/40) - (ext 90/10) 0.010053 0.00349 230 2.880 0.0491
## (ext 60/40) - (acq 75/25) -0.006206 0.00315 271 -1.972 0.3610
## (ext 60/40) - (ext 75/25) -0.003062 0.00315 271 -0.973 0.9261
## (ext 60/40) - (acq 90/10) -0.007117 0.00349 230 -2.039 0.3236
## (ext 60/40) - (ext 90/10) 0.003600 0.00349 230 1.031 0.9070
## (acq 75/25) - (ext 75/25) 0.003144 0.00239 229 1.316 0.7757
## (acq 75/25) - (acq 90/10) -0.000911 0.00310 270 -0.294 0.9997
## (acq 75/25) - (ext 90/10) 0.009806 0.00310 270 3.159 0.0216
## (ext 75/25) - (acq 90/10) -0.004055 0.00310 270 -1.306 0.7814
## (ext 75/25) - (ext 90/10) 0.006662 0.00310 270 2.146 0.2668
## (acq 90/10) - (ext 90/10) 0.010717 0.00345 229 3.103 0.0259
##
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 6 estimateseffectsize::t_to_eta2(t=sumtr$t.ratio, df_error = sumtr$df, type="eta", alternative = "two.sided")## Eta2 (partial) | 95% CI
## -----------------------------
## 0.01 | [0.00, 0.06]
## 2.29e-05 | [0.00, 0.01]
## 4.27e-03 | [0.00, 0.03]
## 1.57e-04 | [0.00, 0.01]
## 0.03 | [0.00, 0.09]
## 0.01 | [0.00, 0.05]
## 3.49e-03 | [0.00, 0.03]
## 0.02 | [0.00, 0.07]
## 4.61e-03 | [0.00, 0.04]
## 7.53e-03 | [0.00, 0.04]
## 3.19e-04 | [0.00, 0.02]
## 0.04 | [0.01, 0.09]
## 6.28e-03 | [0.00, 0.04]
## 0.02 | [0.00, 0.06]
## 0.04 | [0.01, 0.10]# one-way t-tests against true reinforcement rates
df2 <- rev_data %>%
group_by(contingency, tabin, phase_str) %>%
nest() %>%
mutate(tt=map(data,~t.test(.x$err))) %>%
mutate(tidied = map(tt, tidy)) %>%
unnest(tidied, .drop = T)
df2$p.value.fdr = p.adjust(df2$p.value, method = "fdr")
df2$sig <- 0
df2$sig[df2$p.value.fdr<0.05] <- 1
df2## # A tibble: 12 × 15
## # Groups: contingency, tabin, phase_str [12]
## contingency phase_str tabin data tt estimate statistic p.value parameter
## <fct> <fct> <fct> <lis> <lis> <dbl> <dbl> <dbl> <dbl>
## 1 75/25 acq lowA… <tib… <hte… -0.0677 -2.79 7.50e- 3 47
## 2 75/25 ext lowA… <tib… <hte… 0.269 9.42 2.12e-12 47
## 3 75/25 acq hiAnx <tib… <hte… -0.0602 -2.47 1.79e- 2 39
## 4 75/25 ext hiAnx <tib… <hte… 0.214 6.33 1.81e- 7 39
## 5 60/40 acq hiAnx <tib… <hte… -0.0229 -0.662 5.17e- 1 17
## 6 60/40 ext hiAnx <tib… <hte… 0.0936 2.12 4.91e- 2 17
## 7 90/10 acq hiAnx <tib… <hte… -0.103 -2.13 4.80e- 2 17
## 8 90/10 ext hiAnx <tib… <hte… 0.194 3.38 3.59e- 3 17
## 9 60/40 acq lowA… <tib… <hte… -0.0212 -0.502 6.22e- 1 17
## 10 60/40 ext lowA… <tib… <hte… 0.232 5.25 6.47e- 5 17
## 11 90/10 acq lowA… <tib… <hte… -0.125 -3.01 7.47e- 3 18
## 12 90/10 ext lowA… <tib… <hte… 0.351 5.99 1.16e- 5 18
## # … with 6 more variables: conf.low <dbl>, conf.high <dbl>, method <chr>,
## # alternative <chr>, p.value.fdr <dbl>, sig <dbl># binarized anxiety to get means per group for descriptive tables (in supplementary)
m2c <-lmer(prob ~ tabin*contingency*phase_str + (1|id) + (1|study_str), rev_data) ## boundary (singular) fit: see help('isSingular')em2c1 = emmeans(m2c, specs = pairwise ~ phase_str*tabin)## NOTE: Results may be misleading due to involvement in interactionsem2c1$emmeans## phase_str tabin emmean SE df lower.CL upper.CL
## acq lowAnx 0.698 0.0303 6.03 0.624 0.772
## ext lowAnx 0.487 0.0303 6.03 0.412 0.561
## acq hiAnx 0.707 0.0322 6.00 0.628 0.786
## ext hiAnx 0.388 0.0322 6.00 0.309 0.467
##
## Results are averaged over the levels of: contingency
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95em2c2 = emmeans(m2c, specs = pairwise ~ phase_str*contingency*tabin)
em2c2$emmeans## phase_str contingency tabin emmean SE df lower.CL upper.CL
## acq 60/40 lowAnx 0.629 0.0499 28.8 0.527 0.731
## ext 60/40 lowAnx 0.536 0.0499 28.8 0.434 0.638
## acq 75/25 lowAnx 0.712 0.0279 16.6 0.653 0.771
## ext 75/25 lowAnx 0.469 0.0279 16.6 0.410 0.528
## acq 90/10 lowAnx 0.753 0.0489 26.3 0.653 0.854
## ext 90/10 lowAnx 0.455 0.0489 26.3 0.354 0.555
## acq 60/40 hiAnx 0.632 0.0504 27.2 0.528 0.735
## ext 60/40 hiAnx 0.422 0.0504 27.2 0.318 0.525
## acq 75/25 hiAnx 0.715 0.0314 19.3 0.649 0.780
## ext 75/25 hiAnx 0.441 0.0314 19.3 0.376 0.507
## acq 90/10 hiAnx 0.774 0.0504 27.2 0.671 0.878
## ext 90/10 hiAnx 0.301 0.0504 27.2 0.197 0.404
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95em2c3 = emmeans(m2c, specs = pairwise ~ phase_str*contingency)## NOTE: Results may be misleading due to involvement in interactionsem2c3$emmeans## phase_str contingency emmean SE df lower.CL upper.CL
## acq 60/40 0.630 0.0391 7.34 0.539 0.722
## ext 60/40 0.479 0.0391 7.34 0.387 0.571
## acq 75/25 0.713 0.0212 4.80 0.658 0.769
## ext 75/25 0.455 0.0212 4.80 0.400 0.510
## acq 90/10 0.764 0.0388 7.00 0.672 0.856
## ext 90/10 0.378 0.0388 7.00 0.286 0.469
##
## Results are averaged over the levels of: tabin
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95Plots
reversal-locked ratings by state and anxiety
rev_data_cont <- read.csv(here::here("data/reversal_aligned_data.csv"))
rev_data_cont = assign_var_types(rev_data_cont, c("phase_str", "study_str", "id", "contingency", "ta", "half_str"))
#first summarize by participant
rev_data_cont_df = rev_data_cont %>%
group_by(phase_str, contingency, id, trno) %>%
summarise_at("prob", funs(mean),na.rm = TRUE)
df = rev_data_cont_df %>%
group_by(phase_str, trno) %>%
summarise_at("prob", funs(mean,std.error),na.rm = TRUE)
df$lower = df$mean - df$std.error
df$upper = df$mean + df$std.error
g <- ggplot(data = df, aes(y = mean, x = trno, fill=phase_str)) +
geom_segment(x=0.5,y=0.75, xend=15, yend=0.75 , linetype="dashed", color=pal[1], size=0.5 , alpha=0.9) +
geom_segment(x=0.5,y=0.25, xend=15, yend=0.25 , linetype="dashed", color=pal[3], size=0.5 , alpha=0.3) +
geom_segment(x=-5,y=0.75, xend=-0.5, yend=0.75 , linetype="dashed", color=pal[3], size=0.5 , alpha=0.9) +
geom_segment(x=-5,y=0.25, xend=-0.5, yend=0.25 , linetype="dashed", color=pal[1], size=0.5 , alpha=0.3) +
geom_vline(xintercept=0, linetype="dashed") +
geom_line(size=1.5, show.legend = TRUE) +
geom_ribbon(aes(ymin=lower, ymax=upper), na.rm = TRUE,alpha=0.8, show.legend = TRUE) +
geom_rect(xmin=10,xmax=15,ymin=0,ymax=1, size=2, color=pal2[7], fill=rgb(0.9, 0.9, 0.9), alpha=0.01) +
scale_color_manual(values = c(pal[1], pal[3])) +
scale_fill_manual(values = c(pal[1], pal[3]), name = "Switch type", labels = c("Low-to-high (L2H)", "High-to-low (H2L)")) +
scale_x_continuous(breaks = c(-5, -1, 1, 5, 10, 15)) +
ylim(0,1) +
theme_bw() +
labs(y= "Rating", x="Trial (locked to reversal)") +
raincloud_theme +
theme(legend.position = c(0.35,0.9))
g
ggsave(here::here("output/figures/Fig_reversal_1.pdf"), width = 5, height = 5, dpi = 120)reversal-locked ratings by state and anxiety
#first summarize by participant
rev_data_cont_df2 = rev_data_cont %>%
group_by(phase_str, tabin, contingency, id, trno) %>%
summarise_at("prob", funs(mean),na.rm = TRUE)
df = rev_data_cont_df2 %>%
group_by(phase_str,tabin, trno) %>%
summarise_at("prob", funs(mean,std.error),na.rm = TRUE)
df$lower = df$mean - df$std.error
df$upper = df$mean + df$std.error
g <- ggplot(data = df, aes(y = mean, x = trno, fill=interaction(phase_str, tabin))) +
geom_segment(x=0.5,y=0.75, xend=15, yend=0.75 , linetype="dashed", color=pal[1], size=0.5 , alpha=0.9) +
geom_segment(x=0.5,y=0.25, xend=15, yend=0.25 , linetype="dashed", color=pal[3], size=0.5 , alpha=0.3) +
geom_segment(x=-5,y=0.75, xend=-0.5, yend=0.75 , linetype="dashed", color=pal[3], size=0.5 , alpha=0.9) +
geom_segment(x=-5,y=0.25, xend=-0.5, yend=0.25 , linetype="dashed", color=pal[1], size=0.5 , alpha=0.3) +
geom_vline(xintercept=0, linetype="dashed") +
geom_line(size=1.5, show.legend = TRUE) +
geom_ribbon(aes(ymin=lower, ymax=upper), na.rm = TRUE,alpha=0.6,show.legend = TRUE) +
scale_fill_manual(values = c(pal[2], pal[4], pal[1], pal[3]), name = "", labels = c("Low Anx: L2H", "Low Anx: H2L", "High Anx: L2H", "High Anx: H2L")) +
scale_color_manual(values = c(pal[2], pal[4], pal[1], pal[3])) +
ylim(0,1) +
scale_x_continuous(breaks = c(-5, -1, 1, 5, 10), limits=c(-5,10)) +
theme_bw() +
labs(y= "Rating", x="Trial (locked to reversal)") +
raincloud_theme +
theme(legend.position = c(0.75,0.25))
g
ggsave(here::here("output/figures/Fig_reversal_2.pdf"), width = 5, height = 5, dpi = 120)reversal cue by anxiety, session and a state
df = rev_data %>%
group_by(id, phase_str, contingency, tabin) %>%
summarise_at("prob", mean, na.rm = TRUE)
df_levs <- data.frame(contingency = c("60/40", "60/40", "75/25", "75/25", "90/10", "90/10"),
phase_str = c("acq", "ext", "acq", "ext", "acq", "ext"),
level = c(60, 40, 75, 25, 90,10))
g <- ggplot(data = df, aes(y = prob, x = phase_str, fill = interaction( tabin, phase_str))) +
geom_hline(data = df_levs[df_levs$phase_str == "acq",], aes(yintercept = level/100), color=pal[1], linetype="longdash", alpha =0.8) +
geom_hline(data = df_levs[df_levs$phase_str == "ext",], aes(yintercept = level/100), color=pal[3], linetype="longdash", alpha =0.8) +
geom_point(aes(color= interaction( tabin, phase_str)), position = position_jitterdodge( jitter.width = .45, dodge.width = 1), size = 2, alpha = 0.5, show.legend=FALSE) +
geom_boxplot(inherit.aes = TRUE, width = 0.8, lwd=1.2, outlier.shape = NA, alpha = 0.7) +
stat_summary(geom="point", fun="mean", size=4, shape=23, aes(group=interaction(tabin,phase_str)), color="black",
position = position_dodge( 0.8), stroke=1.5, show.legend = FALSE) +
expand_limits(x = 3) +
scale_color_manual(values = c(pal[2], pal[1], pal[4], pal[3])) +
scale_fill_manual(values = c(pal[2], pal[1], pal[4], pal[3]), name = "", labels = c("HS: Low Anxiety", "HS: High Anxiety", "LS: Low Anxiety", "LS: High Anxiety")) +
theme_bw() +
raincloud_theme +
theme(strip.text.x = element_text(size = 14, face = "bold"),
strip.background = element_rect(color="black", fill="#ebebeb", size=0, linetype="solid")) +
theme(legend.position = "top") + #c(0.8, 0.2)
labs(y= "Mean rating", x="State") +
scale_x_discrete(labels=c("acq" = "High", "ext" = "Low")) +
facet_grid(cols=vars(contingency)) +
guides(fill=guide_legend(nrow=2,byrow=TRUE)) +
ylim(0,1)
g
ggsave(here::here("output/figures/Fig_reversal_3.pdf"), width = 6, height = 5, dpi = 120)without session to better reflect the statistical model
df = rev_data %>%
group_by(id, phase_str, contingency, ta_orig_scale) %>%
summarise_at("prob", mean, na.rm = TRUE)
df_levs <- data.frame(contingency = c("75/25", "75/25"),
phase_str = c("acq", "ext", "acq", "ext", "acq", "ext"),
level = c(75, 25))
g <- ggplot(data = df, aes(y = prob, x = ta_orig_scale, fill = interaction(phase_str))) +
geom_hline(data = df_levs[df_levs$phase_str == "acq",], aes(yintercept = level/100), color=pal[1], linetype="longdash", alpha =0.8) +
geom_hline(data = df_levs[df_levs$phase_str == "ext",], aes(yintercept = level/100), color=pal[3], linetype="longdash", alpha =0.8) +
geom_point(aes(color= phase_str), position = position_jitterdodge( jitter.width = .45, dodge.width = 1), size = 2, alpha = 0.5, show.legend=FALSE) +
geom_smooth(aes(color= phase_str, fill=phase_str), method='lm', formula= y~x, alpha=0.3, show.legend = FALSE) +
#geom_boxplot(inherit.aes = TRUE, width = 0.8, lwd=1.2, outlier.shape = NA, alpha = 0.7) +
#stat_summary(geom="point", fun="mean", size=4, shape=23, aes(group=interaction(tabin,phase_str)), color="black",
# position = position_dodge( 0.8), stroke=1.5, show.legend = FALSE) +
expand_limits(x = 3) +
scale_color_manual(values = c(pal[1], pal[3], pal[4], pal[3])) +
scale_fill_manual(values = c(pal[1], pal[3], pal[4], pal[3]), name = "", labels = c("HS: Low Anxiety", "HS: High Anxiety", "LS: Low Anxiety", "LS: High Anxiety")) +
theme_bw() +
raincloud_theme +
theme(strip.text.x = element_text(size = 14, face = "bold"),
strip.background = element_rect(color="black", fill="#ebebeb", size=0, linetype="solid")) +
theme(legend.position = "top") + #c(0.8, 0.2)
labs(y= "Mean rating on trials 10+", x="Trait anxiety") +
#scale_x_discrete(labels=c("acq" = "High", "ext" = "Low")) +
#facet_grid(cols=vars(contingency)) +
guides(fill=guide_legend(nrow=2,byrow=TRUE)) +
ylim(0,1) +
xlim(20,70)
g
ggsave(here::here("output/figures/Fig_reversal_3_collapsed.pdf"), width = 4, height = 5, dpi = 120)Only 75/25 by study
df = rev_data %>%
group_by(id, phase_str, contingency, tabin,ta, study_str) %>%
summarise_at("prob", mean, na.rm = TRUE)
df <- df[df$contingency %in% "75/25",]
m<-lmer(prob ~ phase_str*study_str + (1|id) , df)
anova(m)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## phase_str 2.87820 2.87820 1 85 116.2776 <2e-16 ***
## study_str 0.03662 0.01831 2 85 0.7398 0.4803
## phase_str:study_str 0.05418 0.02709 2 85 1.0945 0.3394
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1joint_tests(m, by = "phase_str")## phase_str = acq:
## model term df1 df2 F.ratio p.value
## study_str 2 158.14 0.351 0.7047
##
## phase_str = ext:
## model term df1 df2 F.ratio p.value
## study_str 2 158.14 1.386 0.2530df_levs <- data.frame(contingency = c("75/25", "75/25", "75/25", "75/25", "75/25", "75/25"),
phase_str = c("acq", "ext", "acq", "ext", "acq", "ext"),
level = c(75, 25, 75, 25, 75, 25))
g <- ggplot(data = df, aes(y = prob, x = phase_str, fill = interaction( tabin, phase_str))) +
geom_hline(data = df_levs[df_levs$phase_str == "acq",], aes(yintercept = level/100), color=pal[1], linetype="longdash", alpha =0.8) +
geom_hline(data = df_levs[df_levs$phase_str == "ext",], aes(yintercept = level/100), color=pal[3], linetype="longdash", alpha =0.8) +
geom_point(aes(color= interaction( tabin, phase_str)), position = position_jitterdodge( jitter.width = .45, dodge.width = 1), size = 2, alpha = 0.5, show.legend=FALSE) +
geom_boxplot(inherit.aes = TRUE, width = 0.8, lwd=1.2, outlier.shape = NA, alpha = 0.7) +
expand_limits(x = 3) +
scale_color_manual(values = c(pal[2], pal[1], pal[4], pal[3])) +
scale_fill_manual(values = c(pal[2], pal[1], pal[4], pal[3]), name = "", labels = c("HS: Low Anxiety", "HS: High Anxiety", "LS: Low Anxiety", "LS: High Anxiety")) +
theme_bw() +
raincloud_theme +
theme(strip.text.x = element_text(size = 14, face = "bold"),
strip.background = element_rect(color="black", fill="#ebebeb", size=0, linetype="solid")) +
theme(legend.position = "top") + #c(0.8, 0.2)
labs(y= "Mean rating", x="State") +
scale_x_discrete(labels=c("acq" = "High", "ext" = "Low")) +
facet_grid(cols=vars(study_str)) +
guides(fill=guide_legend(nrow=2,byrow=TRUE)) +
ylim(0,1)
g
ggsave(here::here("output/figures/SuppFig5_raw.pdf"), width = 6, height = 5, dpi = 120)3. Analysis of slopes and switchpoints
Stats
Slopes by anxiety and state
rev_data <- read.csv(here::here("data/full_REV_data.csv"))
# We'll be fitting the slopes to the first ten trials after reversal
rev_data<-rev_data[which(rev_data$trno %in% seq(10)),]
rev_data = assign_var_types(rev_data, c("phase_str", "trno", "study_str", "id", "contingency", "ta", "half_str"))
# prepare data set
conts <- c("60/40", "75/25", "90/10")
phases <- c("acq", "ext") # acq= high-state; ext = low-state
halves <- c("first", "second")
df <- rev_data[,c("id", "ta")]
df<-df[!duplicated(df$id),]
df["tabin"]<-dicho(df$ta)
df$tabin <- mapvalues(df$tabin,
from = c(0,1),
to = c("lowAnx", "hiAnx"))
df["tabin"]<-to_factor(df["tabin"])
# run lmer for each combination of session, phase and half
dt <- data.frame()
for (c in conts) {
for (ph in phases) {
for (h in halves) {
adata <- rev_data[which((rev_data$phase_str %in% ph) & (rev_data$contingency %in% c) & (rev_data$half_str %in% h)),]
mm = lmer(prob ~ 1 + (trno|id) , adata) # This estimates slope for each participant
rn_a<- data.frame(ranef(mm))
rn_a["id"] <- rn_a$grp
df_a<-data.frame()
df_a <- merge(df, rn_a[rn_a$term=="trno",], by="id")
df_a["phase"] <- ph
df_a["contingency"] <- c
df_a["half"] <- h
dt <- rbind(dt, df_a)
}
}
}
dt$phase_str = to_factor(dt$phase)
dt$contingency = to_factor(dt$contingency)
dt$half_str = to_factor(dt$half)
dt$abscondval <- abs(dt$condval)
# Also save the slope estiamtes for other anlayses
dt$slope <- dt$abscondval
write.csv(dt, here::here("data", "slope_estimates.csv"))
# Ns
N = dt %>%
group_by(contingency) %>%
count()
dt_count <- dt %>%
group_by(contingency) %>%
summarise(count = n_distinct(id))
dt_count## # A tibble: 3 × 2
## contingency count
## <fct> <int>
## 1 60/40 36
## 2 75/25 88
## 3 90/10 37# Test on slopes (absolute slopes)
dfst<-read.csv(here::here("data/model_fit_final_full.csv"))
dfst<-dfst[,c("ids", "study", "contingency")]
dfst<-assign_var_types(dfst, c( "study", "ids"))
dfst$id <- dfst$ids
dfst<-dfst[,c("study_str", "id", "contingency")]
dt <- base::merge(dt, dfst, by=c("id", "contingency"), all.x=TRUE)
# Report mean slopes using non-absolute values
m3b<-lmer(condval ~ ta*contingency*phase_str + (1|id) + (1|study_str), dt)## boundary (singular) fit: see help('isSingular')em3b = emmeans(m3b, specs = pairwise ~ phase_str)## NOTE: Results may be misleading due to involvement in interactionsem3b$emmeans## phase_str emmean SE df lower.CL upper.CL
## acq 0.0244 0.00287 0.75 -0.0578 0.1067
## ext -0.0233 0.00290 0.74 -0.1098 0.0632
##
## Results are averaged over the levels of: contingency
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95jt<-joint_tests(m3b, by = "contingency")
jt## contingency = 60/40:
## model term df1 df2 F.ratio p.value
## ta 1 290.73 0.901 0.3434
## phase_str 1 554.55 33.404 <.0001
## ta:phase_str 1 554.55 4.629 0.0319
##
## contingency = 75/25:
## model term df1 df2 F.ratio p.value
## ta 1 292.03 0.518 0.4725
## phase_str 1 561.13 141.254 <.0001
## ta:phase_str 1 562.70 0.050 0.8234
##
## contingency = 90/10:
## model term df1 df2 F.ratio p.value
## ta 1 290.73 0.005 0.9441
## phase_str 1 554.55 316.719 <.0001
## ta:phase_str 1 554.55 29.456 <.0001# Main statistical model using absolute value of the slope
m3a<-lmer(abscondval ~ ta*contingency*phase_str + (1|id) + (1|study_str), dt)
anova(m3a)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## ta 0.002629 0.0026291 1 87.42 5.8596 0.017561 *
## contingency 0.048342 0.0241710 2 96.11 53.8716 < 2.2e-16 ***
## phase_str 0.000531 0.0005307 1 541.05 1.1828 0.277262
## ta:contingency 0.004787 0.0023936 2 599.71 5.3348 0.005052 **
## ta:phase_str 0.000161 0.0001609 1 542.26 0.3587 0.549489
## contingency:phase_str 0.000383 0.0001914 2 541.29 0.4267 0.652896
## ta:contingency:phase_str 0.000856 0.0004278 2 542.73 0.9535 0.386053
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# statistical assumptions
## homogeneity of variance
dt$residuals <- residuals(m3a)
dt$residuals_sq <- abs(residuals(m3a))^2
hist(dt$residuals_sq)
lev_m <- lm(residuals_sq ~ id, data=dt) #ANOVA of the squared residuals
anova(lev_m)## Analysis of Variance Table
##
## Response: residuals_sq
## Df Sum Sq Mean Sq F value Pr(>F)
## id 88 5.0268e-05 5.7123e-07 1.2109 0.107
## Residuals 544 2.5662e-04 4.7172e-07## normality of residuals
ggplot(data=dt, aes(x=residuals)) +
geom_histogram()## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
car::qqPlot(dt$residuals)
## [1] 244 593# effectsize
effectsize::effectsize(anova(m3a), type="eta", alternative = "two.sided")## # Effect Size for ANOVA (Type III)
##
## Parameter | Eta2 (partial) | 95% CI
## --------------------------------------------------------
## ta | 0.06 | [0.00, 0.18]
## contingency | 0.53 | [0.39, 0.63]
## phase_str | 2.18e-03 | [0.00, 0.02]
## ta:contingency | 0.02 | [0.00, 0.04]
## ta:phase_str | 6.61e-04 | [0.00, 0.01]
## contingency:phase_str | 1.57e-03 | [0.00, 0.01]
## ta:contingency:phase_str | 3.50e-03 | [0.00, 0.02]# # Post-hocs
em3a = emmeans(m3a, specs = pairwise ~ contingency)## NOTE: Results may be misleading due to involvement in interactionsem3a$emmeans## contingency emmean SE df lower.CL upper.CL
## 60/40 0.0200 0.00305 3.72 0.0113 0.0287
## 75/25 0.0252 0.00188 1.85 0.0165 0.0340
## 90/10 0.0445 0.00303 3.57 0.0357 0.0534
##
## Results are averaged over the levels of: phase_str
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95em3a$contrasts## contrast estimate SE df t.ratio p.value
## (60/40) - (75/25) -0.00521 0.00271 88.2 -1.921 0.1387
## (60/40) - (90/10) -0.02452 0.00249 541.2 -9.844 <.0001
## (75/25) - (90/10) -0.01931 0.00269 84.9 -7.191 <.0001
##
## Results are averaged over the levels of: phase_str
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 3 estimatesemstr <- summary(em3a$contrasts)
effectsize::t_to_eta2(t=emstr$t.ratio, df_error = emstr$df, alternative = "two.sided")## Eta2 (partial) | 95% CI
## -----------------------------
## 0.04 | [0.00, 0.15]
## 0.15 | [0.10, 0.21]
## 0.38 | [0.22, 0.51]# # Interactions with trait anxiety
trdf<-emtrends(m3a, pairwise ~ contingency, var = "ta")## NOTE: Results may be misleading due to involvement in interactionstrdf## $emtrends
## contingency ta.trend SE df lower.CL upper.CL
## 60/40 0.000185 0.000208 241 -0.000224 0.000593
## 75/25 0.000110 0.000156 141 -0.000198 0.000417
## 90/10 0.000748 0.000205 233 0.000345 0.001151
##
## Results are averaged over the levels of: phase_str
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## (60/40) - (75/25) 0.000075 0.000213 620 0.353 0.9336
## (60/40) - (90/10) -0.000564 0.000222 542 -2.544 0.0302
## (75/25) - (90/10) -0.000639 0.000209 620 -3.049 0.0067
##
## Results are averaged over the levels of: phase_str
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 3 estimatesemstr <- summary(trdf$contrasts)
effectsize::t_to_eta2(t=emstr$t.ratio, df_error = emstr$df, alternative = "two.sided")## Eta2 (partial) | 95% CI
## -----------------------------
## 2.01e-04 | [0.00, 0.01]
## 0.01 | [0.00, 0.04]
## 0.01 | [0.00, 0.04]jt<-joint_tests(m3a, by = "contingency")
jt## contingency = 60/40:
## model term df1 df2 F.ratio p.value
## ta 1 240.53 0.791 0.3746
## phase_str 1 538.11 0.040 0.8416
## ta:phase_str 1 538.11 0.453 0.5011
##
## contingency = 75/25:
## model term df1 df2 F.ratio p.value
## ta 1 141.34 0.496 0.4823
## phase_str 1 545.42 0.030 0.8622
## ta:phase_str 1 547.94 2.041 0.1536
##
## contingency = 90/10:
## model term df1 df2 F.ratio p.value
## ta 1 232.86 13.385 0.0003
## phase_str 1 538.11 0.691 0.4063
## ta:phase_str 1 538.11 0.379 0.5385jt$F.ratio## [1] 0.791 0.496 13.385 0.040 0.030 0.691 0.453 2.041 0.379effectsize::F_to_eta2(f=jt$F.ratio , df=jt$df1, df_error = jt$df2, alternative = "two.sided")## Eta2 (partial) | 95% CI
## -----------------------------
## 3.28e-03 | [0.00, 0.03]
## 3.50e-03 | [0.00, 0.05]
## 0.05 | [0.01, 0.12]
## 7.43e-05 | [0.00, 0.01]
## 5.50e-05 | [0.00, 0.01]
## 1.28e-03 | [0.00, 0.01]
## 8.41e-04 | [0.00, 0.01]
## 3.71e-03 | [0.00, 0.02]
## 7.04e-04 | [0.00, 0.01]jt2<-joint_tests(m3a, by = c("contingency", "phase_str"))
jt2## contingency = 60/40, phase_str = acq:
## model term df1 df2 F.ratio p.value
## ta 1 435.34 1.244 0.2652
##
## contingency = 75/25, phase_str = acq:
## model term df1 df2 F.ratio p.value
## ta 1 283.04 0.068 0.7950
##
## contingency = 90/10, phase_str = acq:
## model term df1 df2 F.ratio p.value
## ta 1 425.79 6.484 0.0112
##
## contingency = 60/40, phase_str = ext:
## model term df1 df2 F.ratio p.value
## ta 1 435.34 0.091 0.7634
##
## contingency = 75/25, phase_str = ext:
## model term df1 df2 F.ratio p.value
## ta 1 292.30 1.946 0.1641
##
## contingency = 90/10, phase_str = ext:
## model term df1 df2 F.ratio p.value
## ta 1 425.79 10.817 0.0011slope and relative model fit
fit_data <- read.csv(here::here("data/model_fit_final_full.csv"))
fit_data = assign_var_types(fit_data, c("contingency", "ta", "study", "ids"))
fit_data$id <- fit_data$ids
dt$slope <- abs(dt$condval)
dt2 = dt %>%
group_by(id, contingency) %>%
summarise_at(c("slope"), mean, na.rm = TRUE)
write.csv(dt2, here::here("../output", "tables_and_csvs", "est_slopes_1-10.csv"), row.names = FALSE)
df <- join(x=as.data.frame(dt2[,c("contingency", "slope", "id")]), y=as.data.frame(fit_data[, c("contingency", "id", "m1m3_diff", "ta")]), by=c("contingency", "id"))
### correlation between slope and relative model fit
df = df %>%
group_by(id ) %>%
summarise_at(c("slope", "m1m3_diff"), mean, na.rm = TRUE)
df<- as.data.frame(df)
res<-cor.test(df$slope, df$m1m3_diff)
res##
## Pearson's product-moment correlation
##
## data: df$slope and df$m1m3_diff
## t = 3.6114, df = 87, p-value = 0.0005087
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.1652314 0.5295049
## sample estimates:
## cor
## 0.3610637effectsize::t_to_eta2(t=res$statistic, df_error = res$parameter)## Eta2 (partial) | 95% CI
## -----------------------------
## 0.13 | [0.04, 1.00]
##
## - One-sided CIs: upper bound fixed at [1.00].Plots
Slopes by anxiety and state
g <- ggplot(data = dt, aes(y = condval, x = phase, fill = interaction(tabin, phase_str))) +
geom_flat_violin(position = position_nudge(x = .2, y = 0), alpha = .7,lwd=1.2) +
geom_point(aes(color= interaction( tabin, phase_str)), position = position_jitterdodge( jitter.width = .15, dodge.width = 0.2), size = 2, alpha = 0.5, show.legend=FALSE) + geom_boxplot(width = .2, outlier.shape = NA, alpha = 0.7, lwd=1.2, show.legend = FALSE) +
geom_hline(yintercept=c(0), linetype="dashed") +
stat_summary(geom="point", fun="mean", size=3, shape=23, aes(group=interaction(tabin,phase_str)), color="black",
position = position_dodge( 0.2), stroke=1, show.legend = FALSE) +
expand_limits(x = 3) +
scale_color_manual(values = c(pal[2], pal[1], pal[4], pal[3])) +
scale_fill_manual(values = c(pal[2], pal[1], pal[4], pal[3]), name = "Trait anxiety", labels = c("Low", "High")) +
theme_bw() +
raincloud_theme +
theme(strip.text.x = element_text(size = 14, face = "bold"),
strip.background = element_rect(color="black", fill="#ebebeb", size=0, linetype="solid")) +
labs(y= "Slope on trials 1 - 10", x="State") +
scale_x_discrete(labels=c("acq" = "High", "ext" = "Low")) +
scale_y_continuous(labels = function(x) paste0(x*100, "%")) +
facet_grid(cols=vars(contingency))
g
ggsave(here::here("output/figures/Fig_slopes_1.pdf"), width = 10, height = 3, dpi = 120)Correlation of TA and slope in 90/10 session
### assumes previous dt
dt$slope <- abs(dt$condval)
dt2 = dt %>%
group_by(id, contingency) %>%
summarise_at(c("slope"), mean, na.rm = TRUE)
df <- join(x=as.data.frame(dt2[,c("contingency", "slope", "id")]), y=as.data.frame(fit_data[, c("contingency", "id", "m1m3_diff", "ta_orig_scale")]), by=c("contingency", "id"))
df<-df[df$contingency %in% "90/10",]
g <- ggplot(data = df, aes(x = ta_orig_scale, y = slope ), show.legend = FALSE) +
geom_point( ) +
geom_smooth(method='lm', formula= y~x, alpha=0.3, show.legend = FALSE, color="black") +
scale_color_gradient(low=pal[9], high=pal[9]) +
theme_bw() +
raincloud_theme +
stat_cor(method = "pearson", alternative = "two.sided", cor.coef.name = c("r"), label.sep = ", ", label.x.npc = "left", label.y.npc = "top", digits = 3, r.digits = 3, p.digits = 3) +
labs(x= "Trait anxiety", y="Slope on trials 1-10")
g
ggsave(here::here("output/figures/Fig_slopes_2.pdf"), width = 5, height = 4, dpi = 120)Switch steepness and switchpoint
stsw_data <- read.csv(here::here("data/steepness_and_switchpoint_full.csv"))
stsw_data = assign_var_types(stsw_data, c("study", "ta", "phase_str", "half_str", "cont", "tabin", "session"))
stsw_data = stsw_data %>%
group_by(id, cont,phase_str, ta, study_str) %>%
summarise_at(c("log_steepness", "steepness","switchpoint"), mean, na.rm = TRUE)
m <- lmer(log_steepness ~ ta*cont*phase_str + (1|id) + (1|study_str) , stsw_data)
anova(m)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## ta 5.6574 5.6574 1 94.376 7.2956 0.008195 **
## cont 31.0785 15.5392 2 205.307 20.0387 1.122e-08 ***
## phase_str 2.1806 2.1806 1 226.858 2.8120 0.094941 .
## ta:cont 2.3203 1.1602 2 260.232 1.4961 0.225924
## ta:phase_str 0.1585 0.1585 1 226.858 0.2045 0.651579
## cont:phase_str 1.1010 0.5505 2 226.858 0.7099 0.492788
## ta:cont:phase_str 0.1640 0.0820 2 226.858 0.1057 0.899727
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# effectsize
effectsize::effectsize(anova(m), type="eta")## # Effect Size for ANOVA (Type III)
##
## Parameter | Eta2 (partial) | 95% CI
## -------------------------------------------------
## ta | 0.07 | [0.01, 1.00]
## cont | 0.16 | [0.09, 1.00]
## phase_str | 0.01 | [0.00, 1.00]
## ta:cont | 0.01 | [0.00, 1.00]
## ta:phase_str | 9.00e-04 | [0.00, 1.00]
## cont:phase_str | 6.22e-03 | [0.00, 1.00]
## ta:cont:phase_str | 9.31e-04 | [0.00, 1.00]
##
## - One-sided CIs: upper bound fixed at [1.00].col<- "slategray"
### Steepness plot
g <- ggplot(aes(x=cont, y=log_steepness, fill=cont), data=stsw_data) +
geom_flat_violin(position = position_nudge(x = .2, y = 0), lwd=1.2, alpha = .8, show.legend = FALSE) +
geom_point(aes(y = log_steepness, color=cont), position = position_jitter(width = .15), size = 2, alpha = 0.5, show.legend = FALSE) +
geom_boxplot(width = .2, outlier.shape = NA, lwd=1.2,alpha = 0.7, show.legend = FALSE) +
stat_summary(geom="point", fun="mean", size=5, shape=23, aes(group=cont), color="black",
position = position_dodge( 0), stroke=1, show.legend = FALSE) +
expand_limits(x = 3) +
scale_color_manual(values=c(col, col, col)) +
scale_fill_manual(values=c(col, col, col)) +
# coord_flip() +
theme_bw() +
raincloud_theme +
labs(y= "Estimated (log) Steepness", x="Session")
g
ggsave(here::here("output/figures/Fig_steepness.pdf"), width = 5, height = 4, dpi = 120)
### Switchpoitn plot
col <- "slategray"
g <- ggplot(aes(x=cont, y=switchpoint, fill=cont), data=stsw_data) +
geom_flat_violin(position = position_nudge(x = .2, y = 0), lwd=1.2, alpha = .8, show.legend = FALSE) +
geom_point(aes(y = switchpoint, color=cont), position = position_jitter(width = .15), size = 2, alpha = 0.5, show.legend = FALSE) +
geom_boxplot(width = .2, outlier.shape = NA, lwd=1.2, alpha = 0.7, show.legend = FALSE) +
stat_summary(geom="point", fun="mean", size=5, shape=23, aes(group=cont), color="black",
position = position_dodge( 0), stroke=1, show.legend = FALSE) +
expand_limits(x = 3) +
scale_color_manual(values=c(col, col, col)) +
scale_fill_manual(values=c(col, col, col)) +
theme_bw() +
raincloud_theme +
labs(y= "Estimated switch point", x="Session")
g
ggsave(here::here("output/figures/Fig_switchpoint.pdf"), width = 5, height = 4, dpi = 120)4. Analysis of model fit and anxiety
Stats
fit_data <- read.csv(here::here("data/model_fit_final_full.csv"))
fit_data = assign_var_types(fit_data, c("study", "ta","contingency", "ids" ))
# m1m3_diff is relative model fit
m4a<-lmer(data=fit_data, m1m3_diff ~ ta*contingency + (1|ids) + (1|study_str))## boundary (singular) fit: see help('isSingular')anova(m4a)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## ta 659.7 659.72 1 80.511 2.1302 0.148314
## contingency 1413.6 706.82 2 104.817 2.2823 0.107102
## ta:contingency 3217.3 1608.65 2 104.876 5.1942 0.007064 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# statistical assumptions
## homogeneity of variance
fit_data$residuals <- residuals(m4a)
fit_data$residuals_sq <- abs(residuals(m4a))^2
lev_m <- lm(residuals_sq ~ ids, data=fit_data) #ANOVA of the squared residuals
anova(lev_m)## Analysis of Variance Table
##
## Response: residuals_sq
## Df Sum Sq Mean Sq F value Pr(>F)
## ids 88 20178129 229297 1.2707 0.147
## Residuals 72 12992832 180456## normality of residuals
ggplot(data=fit_data, aes(x=residuals)) +
geom_histogram()## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
car::qqPlot(fit_data$residuals)
## [1] 89 68# effectsize
effectsize::effectsize(anova(m4a), type="eta", alternative = "two.sided")## # Effect Size for ANOVA (Type III)
##
## Parameter | Eta2 (partial) | 95% CI
## ----------------------------------------------
## ta | 0.03 | [0.00, 0.13]
## contingency | 0.04 | [0.00, 0.13]
## ta:contingency | 0.09 | [0.01, 0.20]# fit*correlation by session (p-vals uncorrected!)
jt<-joint_tests(m4a, by = "contingency")
jt ## contingency = 60/40:
## model term df1 df2 F.ratio p.value
## ta 1 153.15 0.027 0.8691
##
## contingency = 75/25:
## model term df1 df2 F.ratio p.value
## ta 1 143.10 0.098 0.7542
##
## contingency = 90/10:
## model term df1 df2 F.ratio p.value
## ta 1 153.06 9.611 0.0023effectsize::F_to_eta2(f=jt$F.ratio, df=jt$df1, df_error = jt$df2, type="eta", alternative = "two.sided")## Eta2 (partial) | 95% CI
## -----------------------------
## 1.76e-04 | [0.00, 0.02]
## 6.84e-04 | [0.00, 0.03]
## 0.06 | [0.01, 0.14]# contrasts ta*fit effects between sessions (p-vals adjusted)
em4atr<-emtrends(m4a, pairwise ~ contingency, var = "ta")
em4atr <- summary(em4atr$contrasts)
em4atr## contrast estimate SE df t.ratio p.value
## (60/40) - (75/25) 0.0166 0.338 111.2 0.049 0.9987
## (60/40) - (90/10) -0.9694 0.368 77.5 -2.636 0.0271
## (75/25) - (90/10) -0.9860 0.334 110.5 -2.956 0.0106
##
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 3 estimateseffectsize::t_to_eta2(t=em4atr$t.ratio, df_error = em4atr$df, alternative = "two.sided")## Eta2 (partial) | 95% CI
## -----------------------------
## 2.17e-05 | [0.00, 0.01]
## 0.08 | [0.00, 0.22]
## 0.07 | [0.01, 0.18]Plots
BIC scores (demeaned)
df = fit_data %>%
group_by(ids,contingency,ta) %>%
summarise_at(c("m1_BIC","m3_BIC"), mean, na.rm = TRUE)
# Demean
df <- df %>% gather(model, BIC, m1_BIC:m3_BIC)
dem_df <- df %>%
group_by(contingency) %>%
summarise_at(c("BIC"), mean, na.rm = TRUE)
dem_df$BICmean <- dem_df$BIC
dem_df$BIC <- NaN
df<- merge(x = df, y = dem_df, by = "contingency", all.x = TRUE)
df$BIC_demean <- df$BIC.x - df$BICmean
df_summ = df %>%
group_by(model, contingency) %>%
summarise_each(funs(mean,sd), BIC_demean)
g <- ggplot(df_summ, aes(x=model, y=mean, fill=model)) +
geom_bar(inherit.aes = TRUE, stat = "summary", color="black", show.legend = FALSE, width=0.8, lwd=1) +
scale_fill_manual(values = pal2[c(15,16)]) +
scale_x_discrete(labels=c("1-state", "n-state")) +
labs(y= "BIC (demeaned)", x="") +
facet_grid(~contingency) +
theme_bw() +
raincloud_theme
g## No summary function supplied, defaulting to `mean_se()`
## No summary function supplied, defaulting to `mean_se()`
## No summary function supplied, defaulting to `mean_se()`
ggsave(here::here("output/figures/Fig_model_comparison.pdf"), width = 4, height = 4, dpi = 120)## No summary function supplied, defaulting to `mean_se()`
## No summary function supplied, defaulting to `mean_se()`
## No summary function supplied, defaulting to `mean_se()`Relative model fit and anxiety in 90/10
df <- fit_data[fit_data$contingency %in% c("90/10"),]
g <- ggplot(data = df, aes(x = m1m3_diff, y = ta_orig_scale, color=(m1m3_diff) )) +
geom_point( size = 5, alpha = 1, show.legend = FALSE) +
geom_smooth(method='lm', formula= y~x, alpha=0.3, show.legend = FALSE, color="black", fill="lightgray") +
#scale_color_gradient2(low=pal2[15], midpoint = 0, mid = "grey39", high=pal2[16]) +
scale_color_gradient2(low=pal2[15], midpoint = 30, mid = pal2[16], high=pal2[16]) +
geom_vline(xintercept = 0, linetype = "longdash") +
theme_bw() +
stat_cor(method = "pearson", alternative = "two.sided", cor.coef.name = c("r"), label.sep = ", ", label.x.npc = "left", label.y.npc = "top", digits = 3, r.digits = 3, p.digits = 3) +
raincloud_theme +
labs(y= "Trait anxiety", x="Relative fit")
g
ggsave(here::here("output/figures/Fig_model_fit_and_anxiety.pdf"), width = 7, height = 4.5, dpi = 120)
res<-cor.test(df$ta_orig_scale, df$m1m3_diff)
res##
## Pearson's product-moment correlation
##
## data: df$ta_orig_scale and df$m1m3_diff
## t = 2.5135, df = 35, p-value = 0.01671
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.07673412 0.63464001
## sample estimates:
## cor
## 0.3910308effectsize::t_to_eta2(t=res$statistic, df_error = res$parameter)## Eta2 (partial) | 95% CI
## -----------------------------
## 0.15 | [0.02, 1.00]
##
## - One-sided CIs: upper bound fixed at [1.00].4b. Analysis of winning model parameters
N-state
data <- read.csv(here::here("data/models_no_states_and_switches_data.csv"))
vars <- c("contingency", "ta", "study", "ids")
data <- assign_var_types(data, vars)
df <- data %>%
melt(id.vars = c("ids", "ta", "study_str", "contingency"), measure.vars = c("tau_sh_reversalbsim", "tau_nosh_reversalbsim"), value.name = "tau", variable.name = "outcome" )
df$outcome <- mapvalues(df$outcome,
from = c("tau_sh_reversalbsim", "tau_nosh_reversalbsim"),
to =c("sh", "nosh"))
m5a <- lmer(data=df, tau ~ contingency*outcome*ta + (1|ids) + (1|study_str) )
anova(m5a)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## contingency 0.3056 0.1528 2 92.534 0.5315 0.5895
## outcome 10.6448 10.6448 1 228.265 37.0269 4.873e-09 ***
## ta 0.0615 0.0615 1 82.771 0.2138 0.6450
## contingency:outcome 1.1044 0.5522 2 228.265 1.9208 0.1488
## contingency:ta 0.2530 0.1265 2 273.646 0.4400 0.6445
## outcome:ta 0.0472 0.0472 1 228.265 0.1641 0.6858
## contingency:outcome:ta 0.0249 0.0125 2 228.265 0.0434 0.9576
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 # effectsize
effectsize::effectsize(anova(m5a), type="eta", alternative="two.sided")## # Effect Size for ANOVA (Type III)
##
## Parameter | Eta2 (partial) | 95% CI
## ------------------------------------------------------
## contingency | 0.01 | [0.00, 0.07]
## outcome | 0.14 | [0.07, 0.22]
## ta | 2.58e-03 | [0.00, 0.06]
## contingency:outcome | 0.02 | [0.00, 0.06]
## contingency:ta | 3.21e-03 | [0.00, 0.02]
## outcome:ta | 7.18e-04 | [0.00, 0.02]
## contingency:outcome:ta | 3.80e-04 | [0.00, 0.00]# statistical assumptions
## homogeneity of variance
df$residuals <- residuals(m5a)
df$residuals_sq <- abs(residuals(m5a))^2
lev_m <- lm(residuals_sq ~ ids, data=df) #ANOVA of the squared residuals
anova(lev_m)## Analysis of Variance Table
##
## Response: residuals_sq
## Df Sum Sq Mean Sq F value Pr(>F)
## ids 88 8.6164 0.097914 1.2494 0.09625 .
## Residuals 233 18.2597 0.078368
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## normality of residuals
ggplot(data=df, aes(x=residuals)) +
geom_histogram()## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
car::qqPlot(df$residuals)
## [1] 70 66 memm <- emmeans(m5a, specs = pairwise ~ outcome)## NOTE: Results may be misleading due to involvement in interactions memm$emmeans## outcome emmean SE df lower.CL upper.CL
## sh 1.125 0.0835 2.21 0.797 1.45
## nosh 0.729 0.0835 2.21 0.401 1.06
##
## Results are averaged over the levels of: contingency
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.951-state
data <- read.csv(here::here("data/models_no_states_and_switches_data.csv"))
vars <- c("contingency", "ta", "study", "ids")
data <- assign_var_types(data, vars)
df <- data %>%
melt(id.vars = c("ids", "ta", "study_str", "contingency"), measure.vars = c("tau_sh_reversallbetav5", "tau_nosh_reversallbetav5"), value.name = "tau", variable.name = "outcome" )
df$outcome <- mapvalues(df$outcome,
from = c("tau_sh_reversalbssm", "tau_nosh_reversalbssm"),
to =c("sh", "nosh"))## The following `from` values were not present in `x`: tau_sh_reversalbssm, tau_nosh_reversalbssm m <- lmer(data=df, tau ~ contingency*outcome*ta + (1|ids) + (1|study_str) )
anova(m)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## contingency 1.3266 0.6633 2 164.13 2.5281 0.08291 .
## outcome 6.5164 6.5164 1 227.86 24.8369 1.234e-06 ***
## ta 0.1749 0.1749 1 91.99 0.6665 0.41639
## contingency:outcome 0.6815 0.3408 2 227.86 1.2988 0.27488
## contingency:ta 0.3015 0.1508 2 265.93 0.5746 0.56363
## outcome:ta 0.0854 0.0854 1 227.86 0.3256 0.56882
## contingency:outcome:ta 0.0698 0.0349 2 227.86 0.1330 0.87552
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 # effectsize
effectsize::effectsize(anova(m), type="eta", alternative="two.sided")## # Effect Size for ANOVA (Type III)
##
## Parameter | Eta2 (partial) | 95% CI
## ------------------------------------------------------
## contingency | 0.03 | [0.00, 0.09]
## outcome | 0.10 | [0.04, 0.18]
## ta | 7.19e-03 | [0.00, 0.08]
## contingency:outcome | 0.01 | [0.00, 0.05]
## contingency:ta | 4.30e-03 | [0.00, 0.03]
## outcome:ta | 1.43e-03 | [0.00, 0.03]
## contingency:outcome:ta | 1.17e-03 | [0.00, 0.01] memm <- emmeans(m, specs = pairwise ~ outcome)## NOTE: Results may be misleading due to involvement in interactions memm$emmeans## outcome emmean SE df lower.CL upper.CL
## tau_sh_reversallbetav5 1.166 0.0915 2.44 0.832 1.50
## tau_nosh_reversallbetav5 0.856 0.0915 2.44 0.522 1.19
##
## Results are averaged over the levels of: contingency
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95 memm <- emmeans(m, specs = pairwise ~ outcome*contingency)## NOTE: Results may be misleading due to involvement in interactions memm$emmeans## outcome contingency emmean SE df lower.CL upper.CL
## tau_sh_reversallbetav5 60/40 1.037 0.1255 7.60 0.744 1.33
## tau_nosh_reversallbetav5 60/40 0.857 0.1255 7.60 0.565 1.15
## tau_sh_reversallbetav5 75/25 1.172 0.0866 3.01 0.897 1.45
## tau_nosh_reversallbetav5 75/25 0.761 0.0866 3.01 0.486 1.04
## tau_sh_reversallbetav5 90/10 1.289 0.1244 7.31 0.997 1.58
## tau_nosh_reversallbetav5 90/10 0.950 0.1244 7.31 0.658 1.24
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.955. Analysis of model-free learning rates from oddball and meaningful events
Stats
Oddball/meaninful by anx and outcome
# Load calcualted model-free learning rates
mflr_data <- read.csv(here::here("data/oddball_data_n5.csv"))
mflr_data = assign_var_types(mflr_data, c("study", "lrtype", "ta","contingency", "ids" ))
df <- melt(mflr_data, measure.vars=c("mflr_sh","mflr_nosh"),
value.name = "mflr", variable.name="outcome_str")
df = assign_var_types(df, c("outcome_str" ))
## LOG learningr ates for analysis
df$log_mflr <- log(df$mflr)
df_stats<-df[!is.na(df$log_mflr) & !is.infinite(df$log_mflr),]
m6a<-lmer(log_mflr ~ ta*contingency*lrtype*outcome_str + (1|ids) + (1|study_str), df_stats)## boundary (singular) fit: see help('isSingular')anova(m6a)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value
## ta 0.0001 0.0001 1 97.48 0.0004
## contingency 1.7758 0.8879 2 587.42 2.2837
## lrtype 7.0159 7.0159 1 533.22 18.0453
## outcome_str 16.0828 16.0828 1 533.30 41.3657
## ta:contingency 0.8964 0.4482 2 589.62 1.1527
## ta:lrtype 3.2520 3.2520 1 533.52 8.3643
## contingency:lrtype 1.9219 0.9609 2 533.20 2.4716
## ta:outcome_str 0.9197 0.9197 1 533.67 2.3654
## contingency:outcome_str 1.7216 0.8608 2 533.21 2.2140
## lrtype:outcome_str 0.0794 0.0794 1 533.30 0.2043
## ta:contingency:lrtype 1.6488 0.8244 2 533.30 2.1204
## ta:contingency:outcome_str 1.5256 0.7628 2 533.25 1.9620
## ta:lrtype:outcome_str 0.3092 0.3092 1 533.67 0.7954
## contingency:lrtype:outcome_str 0.0866 0.0433 2 533.21 0.1114
## ta:contingency:lrtype:outcome_str 0.1033 0.0517 2 533.25 0.1328
## Pr(>F)
## ta 0.984937
## contingency 0.102810
## lrtype 2.546e-05 ***
## outcome_str 2.810e-10 ***
## ta:contingency 0.316483
## ta:lrtype 0.003983 **
## contingency:lrtype 0.085416 .
## ta:outcome_str 0.124645
## contingency:outcome_str 0.110262
## lrtype:outcome_str 0.651486
## ta:contingency:lrtype 0.120988
## ta:contingency:outcome_str 0.141589
## ta:lrtype:outcome_str 0.372888
## contingency:lrtype:outcome_str 0.894576
## ta:contingency:lrtype:outcome_str 0.875626
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# effectsize
effectsize::effectsize(anova(m6a), type="eta", alternative="two.sided")## # Effect Size for ANOVA (Type III)
##
## Parameter | Eta2 (partial) | 95% CI
## -----------------------------------------------------------------
## ta | 3.68e-06 | [0.00, 0.00]
## contingency | 7.72e-03 | [0.00, 0.03]
## lrtype | 0.03 | [0.01, 0.07]
## outcome_str | 0.07 | [0.04, 0.12]
## ta:contingency | 3.89e-03 | [0.00, 0.02]
## ta:lrtype | 0.02 | [0.00, 0.04]
## contingency:lrtype | 9.19e-03 | [0.00, 0.03]
## ta:outcome_str | 4.41e-03 | [0.00, 0.02]
## contingency:outcome_str | 8.24e-03 | [0.00, 0.03]
## lrtype:outcome_str | 3.83e-04 | [0.00, 0.01]
## ta:contingency:lrtype | 7.89e-03 | [0.00, 0.03]
## ta:contingency:outcome_str | 7.30e-03 | [0.00, 0.03]
## ta:lrtype:outcome_str | 1.49e-03 | [0.00, 0.01]
## contingency:lrtype:outcome_str | 4.18e-04 | [0.00, 0.01]
## ta:contingency:lrtype:outcome_str | 4.98e-04 | [0.00, 0.01]# statistical assumptions
## homogeneity of variance
df_stats$residuals <- residuals(m6a)
df_stats$residuals_sq <- abs(residuals(m6a))^2
lev_m <- lm(residuals_sq ~ ids, data=df_stats) #ANOVA of the squared residuals
anova(lev_m)## Analysis of Variance Table
##
## Response: residuals_sq
## Df Sum Sq Mean Sq F value Pr(>F)
## ids 88 45.37 0.51556 1.4306 0.009747 **
## Residuals 544 196.05 0.36038
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1levRes = car::leveneTest(residuals_sq ~ ids, data=df_stats)
levRes## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 88 1.4094 0.01266 *
## 544
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## normality of residuals
ggplot(data=df_stats, aes(x=residuals)) +
geom_histogram()## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
car::qqPlot(df_stats$residuals)
## [1] 236 80# post hocs
em5a1 = emmeans(m6a, specs = pairwise ~ contingency:outcome_str)## NOTE: Results may be misleading due to involvement in interactionsem5a1$emmeans## contingency outcome_str emmean SE df lower.CL upper.CL
## 60/40 mflr_sh -1.49 0.1035 7.45 -1.73 -1.25
## 75/25 mflr_sh -1.40 0.0603 3.92 -1.57 -1.23
## 90/10 mflr_sh -1.52 0.1033 7.38 -1.76 -1.28
## 60/40 mflr_nosh -1.67 0.1035 7.45 -1.91 -1.43
## 75/25 mflr_nosh -1.84 0.0606 3.97 -2.01 -1.67
## 90/10 mflr_nosh -1.94 0.1028 7.23 -2.18 -1.70
##
## Results are averaged over the levels of: lrtype
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95em5a1$contrasts## contrast estimate SE df t.ratio p.value
## (60/40 mflr_sh) - (75/25 mflr_sh) -0.0876 0.1037 130 -0.845 0.9585
## (60/40 mflr_sh) - (90/10 mflr_sh) 0.0291 0.1043 527 0.279 0.9998
## (60/40 mflr_sh) - (60/40 mflr_nosh) 0.1829 0.1042 525 1.755 0.4963
## (60/40 mflr_sh) - (75/25 mflr_nosh) 0.3473 0.1036 135 3.353 0.0130
## (60/40 mflr_sh) - (90/10 mflr_nosh) 0.4530 0.1039 527 4.360 0.0002
## (75/25 mflr_sh) - (90/10 mflr_sh) 0.1166 0.1034 130 1.128 0.8691
## (75/25 mflr_sh) - (60/40 mflr_nosh) 0.2704 0.1037 130 2.609 0.1024
## (75/25 mflr_sh) - (75/25 mflr_nosh) 0.4348 0.0674 526 6.451 <.0001
## (75/25 mflr_sh) - (90/10 mflr_nosh) 0.5406 0.1030 128 5.250 <.0001
## (90/10 mflr_sh) - (60/40 mflr_nosh) 0.1538 0.1043 527 1.474 0.6808
## (90/10 mflr_sh) - (75/25 mflr_nosh) 0.3182 0.1033 134 3.080 0.0296
## (90/10 mflr_sh) - (90/10 mflr_nosh) 0.4240 0.1037 526 4.087 0.0007
## (60/40 mflr_nosh) - (75/25 mflr_nosh) 0.1644 0.1036 135 1.588 0.6081
## (60/40 mflr_nosh) - (90/10 mflr_nosh) 0.2702 0.1039 527 2.600 0.0990
## (75/25 mflr_nosh) - (90/10 mflr_nosh) 0.1058 0.1028 132 1.028 0.9076
##
## Results are averaged over the levels of: lrtype
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 6 estimatesem5a2 = emmeans(m6a, specs = pairwise ~ outcome_str)## NOTE: Results may be misleading due to involvement in interactionsem5a2$emmeans## outcome_str emmean SE df lower.CL upper.CL
## mflr_sh -1.47 0.0689 1.78 -1.80 -1.14
## mflr_nosh -1.82 0.0690 1.78 -2.15 -1.48
##
## Results are averaged over the levels of: contingency, lrtype
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95em5a2$contrasts## contrast estimate SE df t.ratio p.value
## mflr_sh - mflr_nosh 0.347 0.0539 526 6.439 <.0001
##
## Results are averaged over the levels of: contingency, lrtype
## Degrees-of-freedom method: kenward-roger# relationship with anxiety
emtrm6<-emtrends(m6a, pairwise ~ lrtype, var = "ta")## NOTE: Results may be misleading due to involvement in interactionsemtrm6 <- summary(emtrm6$contrasts)
emtrm6## contrast estimate SE df t.ratio p.value
## common - odd 0.0141 0.00487 526 2.892 0.0040
##
## Results are averaged over the levels of: contingency, outcome_str
## Degrees-of-freedom method: kenward-rogereffectsize::t_to_eta2(t=emtrm6$t.ratio, df_error = emtrm6$df, alternative = "two.sided")## Eta2 (partial) | 95% CI
## -----------------------------
## 0.02 | [0.00, 0.04]Oddball/meaningful by best model per participant (LOG)
data <- read.csv(here::here("data/oddball_data_n5.csv"))
data$mflr_overall_log <- log(data$mflr_overall)
data$lrtype = to_factor(data$lrtype)
data$ids = to_factor(data$ids)
data$ta <- demean(data$ta)
data["tabin"]<-dicho(data$ta)
data$tabin <- mapvalues(data$tabin,
from = c(0,1),
to = c("lowAnx", "hiAnx"))
odd<-data[data$lrtype %in% "odd",]
common<-data[data$lrtype %in% "common",]
data <- merge(odd, common, by="specID", all.x = TRUE)
data$mflr_diff <- (data$mflr_overall.y) - (data$mflr_overall.x)
data$mflr_diff_log <- (data$mflr_overall_log.y) - (data$mflr_overall_log.x)
data<-data[!is.na(data$mflr_diff_log) & !is.infinite(data$mflr_diff_log),]
data <- data[, -grep(".y", colnames(data))]
colnames(data)<-gsub(".x","",colnames(data))
data <- data %>%
group_by(ids) %>%
summarise_at(c("mflr_diff_log", "mflr_diff"), mean, na.rm = TRUE)
fitd<-read.csv(here::here("data/model_fit_final_full.csv"))
fitd<-assign_var_types(fitd, c("contingency", "ta", "study", "ids"))
fitd <- fitd %>%
group_by(ids, ta, study_str) %>%
summarise_at(c("m1_BIC", "m3_BIC", "m1m3_diff"), mean, na.rm = TRUE)
fitd["best_model"] <- apply(fitd[,c("m1_BIC", "m3_BIC")], 1, which.min)
fitd$best_model <- mapvalues(fitd$best_model,
from = c(1,2),
to =c("1-state", "n-state"))
fitd <- fitd[,c("ids", "best_model", "study_str", "m1m3_diff")]
data <-merge(data, fitd, by="ids")
## analysis
data$best_model <- as.factor(data$best_model)
tt = t.test(mflr_diff_log~best_model, alternative = "two.sided", var.equal = FALSE, data=data)
tt##
## Welch Two Sample t-test
##
## data: mflr_diff_log by best_model
## t = -2.107, df = 68.788, p-value = 0.03877
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.321982745 -0.008787365
## sample estimates:
## mean in group 1-state mean in group n-state
## 0.1295605 0.2949456effectsize::t_to_eta2(t=tt$statistic, df_error = tt$parameter, alternative = "two.sided")## Eta2 (partial) | 95% CI
## -----------------------------
## 0.06 | [0.00, 0.20] data %>%
group_by(best_model) %>%
summarise_at(c("mflr_diff_log", "mflr_diff"), mean, na.rm = TRUE)## # A tibble: 2 × 3
## best_model mflr_diff_log mflr_diff
## <fct> <dbl> <dbl>
## 1 1-state 0.130 0.0212
## 2 n-state 0.295 0.0586cor.test(data$mflr_diff_log, data$m1m3_diff, method = "spearman")##
## Spearman's rank correlation rho
##
## data: data$mflr_diff_log and data$m1m3_diff
## S = 90596, p-value = 0.03121
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.228839confintr::ci_cor(data$mflr_diff_log, data$m1m3_diff, method = "spearman",type="bootstrap")##
## Two-sided 95% bootstrap confidence interval for the true Spearman
## correlation coefficient based on 9999 bootstrap replications and the
## bca method
##
## Sample estimate: 0.228839
## Confidence interval:
## 2.5% 97.5%
## 0.005350783 0.428545212odd_model_data = data
write.csv(odd_model_data, here::here("data", "oddball_by_model.csv"))Oddball/meaningful by anxiety (LOG)
df2 <- df %>%
group_by(tabin, lrtype, ids, ta_orig_scale) %>%
summarise_at(c("log_mflr"), mean, na.rm = TRUE)
dfnew <- df2 %>% tidyr::pivot_wider(
id_cols=c("ids", , "tabin", "ta_orig_scale"),
names_from=lrtype,
values_from=log_mflr
)
dfnew$mflr_diff <- dfnew$common - dfnew$odd
dfnew<-dfnew[!is.na(dfnew$mflr_diff) & !is.infinite(dfnew$mflr_diff),]
# somewhat artificially also do t-test
tt = t.test(mflr_diff~tabin, alternative = "two.sided", var.equal = FALSE, data=dfnew)
tt##
## Welch Two Sample t-test
##
## data: mflr_diff by tabin
## t = -2.2445, df = 77.049, p-value = 0.02767
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.34770972 -0.02078966
## sample estimates:
## mean in group lowAnx mean in group hiAnx
## 0.1691060 0.3533557effectsize::t_to_eta2(t=tt$statistic, df_error = tt$parameter, alternative = "two.sided")## Eta2 (partial) | 95% CI
## -----------------------------
## 0.06 | [0.00, 0.19]dfnew %>%
group_by(tabin) %>%
summarise_at(c("mflr_diff"), mean, na.rm = TRUE)## # A tibble: 2 × 2
## tabin mflr_diff
## <fct> <dbl>
## 1 lowAnx 0.169
## 2 hiAnx 0.353cor.test(dfnew$mflr_diff, dfnew$ta_orig_scale, method = "spearman")##
## Spearman's rank correlation rho
##
## data: dfnew$mflr_diff and dfnew$ta_orig_scale
## S = 81629, p-value = 0.03324
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.2298771confintr::ci_cor(dfnew$mflr_diff, dfnew$ta_orig_scale, method = "spearman",type="bootstrap")##
## Two-sided 95% bootstrap confidence interval for the true Spearman
## correlation coefficient based on 9999 bootstrap replications and the
## bca method
##
## Sample estimate: 0.2298771
## Confidence interval:
## 2.5% 97.5%
## 0.03696545 0.41365519Plots
Oddball/Meaninful by anxiety
df2 <- df %>%
group_by(tabin, lrtype, ids, ta_orig_scale) %>%
summarise_at(c("mflr"), mean, na.rm = TRUE)
dfnew <- df2 %>% tidyr::pivot_wider(
id_cols=c("ids", , "tabin", "ta_orig_scale"),
names_from=lrtype,
values_from=mflr
)
dfnew$mflr_diff <- dfnew$common - dfnew$odd
g <- ggplot(data = dfnew, aes(y = mflr_diff, x = tabin, fill = tabin)) +
geom_flat_violin(position = position_nudge(x = .2, y = 0), alpha = .8, lwd=1.2, show.legend = TRUE) +
geom_point(aes(color= tabin), position = position_jitterdodge( jitter.width = .15, dodge.width = 0.25), size = 2, alpha = 0.5, show.legend=FALSE) +
geom_boxplot(width = .25, outlier.shape = NA, alpha = 0.7, lwd=1.2, show.legend = FALSE) +
stat_summary(geom="point", fun="mean", size=5, shape=23, aes(group=tabin), color="black",
position = position_dodge( 0.25), stroke=1, show.legend = FALSE) +
expand_limits(x = 3) +
#scale_color_manual(values = pal3[c(2,4)] , name = "", labels = c("Low Anx", "High Anx")) +
#scale_fill_manual(values = pal3[c(2,4)], name = "", labels = c("Low Anx", "High Anx")) +
scale_color_manual(values = c("olivedrab2", "darkgreen") , name = "", labels = c("Low Anx", "High Anx")) +
scale_fill_manual(values = c("olivedrab2", "darkgreen"), name = "", labels = c("Low Anx", "High Anx")) +
theme_bw() +
raincloud_theme +
theme(strip.text.x = element_text(size = 14, face = "bold"),
strip.background = element_rect(color="black", fill="#ebebeb", size=0, linetype="solid")) +
theme(legend.position = "top") +
labs(y= "Meaningful-oddball learning difference", x="Trait anxiety") +
scale_x_discrete(labels=c("Low", "High")) +
guides(fill=guide_legend(nrow=4,byrow=TRUE)) +
ylim(-0.1, 0.35)
g
ggsave(here::here("output/figures/Fig_oddball_anxiety.pdf"), width = 4, height = 5, dpi = 120)Also plot continously
df2 <- df %>%
group_by(tabin, lrtype, ids, ta_orig_scale) %>%
summarise_at(c("log_mflr"), mean, na.rm = TRUE)
dfnew <- df2 %>% tidyr::pivot_wider(
id_cols=c("ids", , "tabin", "ta_orig_scale"),
names_from=lrtype,
values_from=log_mflr
)
dfnew$mflr_diff <- dfnew$common - dfnew$odd
dfnew<-dfnew[!is.na(dfnew$mflr_diff) & !is.infinite(dfnew$mflr_diff),]
g <- ggplot(data = dfnew, aes(x = ta_orig_scale, y = mflr_diff, color=ta_orig_scale ), show.legend = FALSE) +
geom_point(size=3 ) +
geom_smooth(method='lm', formula= y~x, alpha=0.3, show.legend = FALSE, color="black") +
scale_color_gradient2(low="olivedrab2", high="darkgreen", midpoint = 40, mid = "olivedrab") +
#scale_color_gradient2(low=pal2[15], midpoint = 0, mid = "grey39", high=pal2[16]) +
theme_bw() +
raincloud_theme +
stat_cor(method = "spearman", alternative = "two.sided", cor.coef.name = c("r"), label.sep = ", ", label.x.npc = "left", label.y.npc = "top", digits = 3, r.digits = 3, p.digits = 3) +
labs(x= "Trait anxiety", y="Meaningful-oddball learning difference")
g
ggsave(here::here("output/figures/Fig_oddball_anxiety_continuous.pdf"), width = 6.5, height = 5, dpi = 120)df2 <- df %>%
group_by(contingency, tabin, lrtype, ids) %>%
summarise_at(c("log_mflr"), mean, na.rm = TRUE)
g <- ggplot(data = df2, aes(y = log_mflr, x = tabin, fill = interaction(tabin, lrtype))) +
geom_flat_violin(position = position_nudge(x = .2, y = 0), alpha = .8, lwd=1.2, show.legend = TRUE) +
geom_point(aes(color= interaction( tabin, lrtype)), position = position_jitterdodge( jitter.width = .15, dodge.width = 0.25), size = 2, alpha = 0.5, show.legend=FALSE) +
geom_boxplot(width = .2, outlier.shape = NA, alpha = 0.7, lwd=1.2, show.legend = FALSE) +
stat_summary(geom="point", fun="mean", size=5, shape=23, aes(group= interaction(tabin, lrtype)), color="black",
position = position_dodge( 0.25), stroke=1, show.legend = FALSE) +
expand_limits(x = 3) +
scale_color_manual(values = pal[c(2,4,1,3,2,4)] ) +
scale_fill_manual(values = pal[c(2,4,1,3,2,4)]) +
theme_bw() +
raincloud_theme +
theme(legend.position = "right") +
labs(y= "Learning rate", x="") +
scale_x_discrete(labels=c("common" = "Meaningful", "odd" = "Oddball")) +
guides(fill=guide_legend(nrow=4,byrow=TRUE))
g
Oddball/Meaninful by best fitted model
## Plot the difference as well
g <- ggplot(data = odd_model_data, aes(y = mflr_diff, x = best_model, fill = best_model)) +
geom_flat_violin(position = position_nudge(x = .2, y = 0), alpha = .8, lwd=1.2) +
geom_point(aes(color= best_model), position = position_jitterdodge( jitter.width = .45, dodge.width = 1), size = 2, alpha = 0.5, show.legend=FALSE) +
geom_boxplot(inherit.aes = TRUE, width = 0.25, lwd=1.2, outlier.shape = NA, alpha = 0.7, show.legend = FALSE) +
stat_summary(geom="point", fun="mean", size=5, shape=23, aes(group=best_model), color="black",
position = position_dodge( 1), stroke=1, show.legend = FALSE) +
expand_limits(x = 3) +
scale_color_manual(values = c(pal2[15], pal2[16], pal[4], pal[3])) +
scale_fill_manual(values = c(pal2[15], pal2[16], pal[4], pal[3]), name = "", labels = c("1-state", "n-state")) +
theme_bw() +
raincloud_theme +
theme(strip.text.x = element_text(size = 14, face = "bold"),
strip.background = element_rect(color="black", fill="#ebebeb", size=0, linetype="solid")) +
theme(legend.position = "top") + #c(0.8, 0.2)
labs(y= "Meaningful-oddball learning difference", x="Best model") +
guides(fill=guide_legend(nrow=1,byrow=TRUE)) +
ylim(-0.1, 0.35)
g
ggsave(here::here("output/figures/Fig_oddball_model.pdf"), width = 4, height = 5, dpi = 120)g <- ggplot(data = odd_model_data, aes(x = m1m3_diff, y = mflr_diff_log, color=m1m3_diff ), show.legend = FALSE) +
geom_point(size=3 ) +
geom_smooth(method='lm', formula= y~x, alpha=0.3, show.legend = FALSE, color="black") +
#scale_color_gradient2(low=pal[3], high=pal[4], midpoint = -5, mid = pal[3]) +
scale_color_gradient2(low=pal2[15], midpoint = 30, mid = pal2[16], high=pal2[16]) +
#scale_color_gradient2(low=pal2[15], midpoint = 0, mid = "grey39", high=pal2[16]) +
geom_vline(xintercept = 0, linetype = "longdash") +
theme_bw() +
raincloud_theme +
stat_cor(method = "spearman", alternative = "two.sided", cor.coef.name = c("r"), label.sep = ", ", label.x.npc = "left", label.y.npc = "top", digits = 3, r.digits = 3, p.digits = 3) +
labs(x= "Model difference", y="Meaningful-oddball learning difference")
g
ggsave(here::here("output/figures/Fig_oddball_model_continuous.pdf"), width = 6.5, height = 5, dpi = 120)